diff --git a/GeoModelExamples/SiliconSystemExample/src/SiliconSystemPlugin.cxx b/GeoModelExamples/SiliconSystemExample/src/SiliconSystemPlugin.cxx
index abfeb3933691a5668750ebc0a1e55d264796beaf..9b56e4c9f612da3238cd2a331bbd479985b78d25 100644
--- a/GeoModelExamples/SiliconSystemExample/src/SiliconSystemPlugin.cxx
+++ b/GeoModelExamples/SiliconSystemExample/src/SiliconSystemPlugin.cxx
@@ -119,18 +119,35 @@ void SiliconSystemPlugin::create(GeoVPhysVol *world, bool /*publish*/) {
   GeoLogVol  *boxLog = new GeoLogVol("BoxLog",boxShape,Air);
   GeoPhysVol *worldBOX = new GeoPhysVol(boxLog);
 
-  for (int k=0;k<3;k++) {
-    for (int j=0;j<3;j++)  {
-      GeoBox *box = new GeoBox (0.1,(j+1)*1.7,20);
-      GeoLogVol  *boxLog=new GeoLogVol("SiDet", box, Aluminium);
+  GeoBox *box1 = new GeoBox (0.1,1*1.7,20);
+  GeoLogVol  *boxLog1 = new GeoLogVol("SiDet", box1, Aluminium);
+
+  GeoBox *box2 = new GeoBox (0.1,2*1.7,20);
+  GeoLogVol  *boxLog2 = new GeoLogVol("SiDet", box2, Aluminium);
 
-      GeoRectSurface* rectSurface = new GeoRectSurface((j+1)*1.9, 20.2);
-      GeoVSurface* surf = new GeoVSurface(rectSurface);
+  GeoBox *box3 = new GeoBox (0.1,3*1.7,20);
+  GeoLogVol  *boxLog3 = new GeoLogVol("SiDet", box3, Aluminium);    
 
+  GeoRectSurface* rectSurface1 = new GeoRectSurface(1*1.9, 20.2);
+  GeoRectSurface* rectSurface2 = new GeoRectSurface(2*1.9, 20.2);
+  GeoRectSurface* rectSurface3 = new GeoRectSurface(3*1.9, 20.2);
+
+  for (int k=0;k<3;k++) {
+    for (int j=0;j<3;j++)  {
+      GeoLogVol* boxLog;
+      if(j==0) boxLog = boxLog1;
+      else if(j==1) boxLog = boxLog2;
+      else boxLog = boxLog3;
+
+      GeoRectSurface* rectSurface;
+      if(j==0) rectSurface = rectSurface1;
+      else if(j==1) rectSurface = rectSurface2;
+      else rectSurface = rectSurface3;
+      
       for (int i=0;i<16;i++) {
         double theta = i/16.0*2*M_PI;
         GeoFullPhysVol *boxPhys=new GeoFullPhysVol(boxLog);
-
+        GeoVSurface* surf = new GeoVSurface(rectSurface);
         // Initial transform is very tricky. Because the Virtual surface is always facing to Z axis initially.
         GeoAlignableTransform* xf0 = new GeoAlignableTransform(GeoTrf::RotateZ3D(theta)*GeoTrf::TranslateX3D((j+1)*8.0+0.4*(i%2))*GeoTrf::TranslateX3D(0.5)*GeoTrf::TranslateZ3D((k-1)*44.0)*GeoTrf::RotateX3D(M_PI/2.0)*GeoTrf::RotateY3D(M_PI/2.0));
         GeoAlignableTransform *xf = new GeoAlignableTransform(GeoTrf::RotateZ3D(theta)*GeoTrf::TranslateX3D((j+1)*8.0+0.4*(i%2)));
@@ -163,14 +180,14 @@ void SiliconSystemPlugin::create(GeoVPhysVol *world, bool /*publish*/) {
     }
     
   }
+  double L=10.0;
+  GeoTrapezoidSurface* trapezoid = new GeoTrapezoidSurface(1.2, 5.2, L); 
+  GeoTrd *trd=new GeoTrd(.2, .2, 1, 5, L);
+  GeoLogVol *trdLog=new GeoLogVol("SiDetEnd", trd,Aluminium);
+
   for (int j=0;j<16;j++) {
     for (int i=-1;i<2;i+=2)  {     
-
-      double L=10.0;
-      GeoTrapezoidSurface* trapezoid = new GeoTrapezoidSurface(1.2, 5.2, L);
       GeoVSurface* surf = new GeoVSurface(trapezoid);
-      GeoTrd *trd=new GeoTrd(.2, .2, 1, 5, L);
-      GeoLogVol *trdLog=new GeoLogVol("SiDetEnd", trd,Aluminium);
       GeoFullPhysVol *trdPhys=new GeoFullPhysVol(trdLog);
 
       GeoTransform *sf0=new GeoTransform(GeoTrf::RotateY3D(j*2*M_PI/16.0)*GeoTrf::RotateX3D(M_PI/2.0)*GeoTrf::TranslateY3D(1.5*L));
diff --git a/documentation/docs/components/kernel/reference/RCBase/GeoVSurfaceShape/Annulus.png b/documentation/docs/components/kernel/reference/RCBase/GeoVSurfaceShape/Annulus.png
new file mode 100644
index 0000000000000000000000000000000000000000..d21fbb3db11f1afc5ff4518480f23d79dcd52931
Binary files /dev/null and b/documentation/docs/components/kernel/reference/RCBase/GeoVSurfaceShape/Annulus.png differ
diff --git a/documentation/docs/components/kernel/reference/RCBase/GeoVSurfaceShape/AnnulusDemo.png b/documentation/docs/components/kernel/reference/RCBase/GeoVSurfaceShape/AnnulusDemo.png
new file mode 100644
index 0000000000000000000000000000000000000000..9a9826b16a555bb28c6d1383df83cd87604759e6
Binary files /dev/null and b/documentation/docs/components/kernel/reference/RCBase/GeoVSurfaceShape/AnnulusDemo.png differ
diff --git a/documentation/docs/components/kernel/reference/RCBase/GeoVSurfaceShape/AnnulusDemo2.png b/documentation/docs/components/kernel/reference/RCBase/GeoVSurfaceShape/AnnulusDemo2.png
new file mode 100644
index 0000000000000000000000000000000000000000..e26e38f03c1546a914f4244162b2aa2c161ec278
Binary files /dev/null and b/documentation/docs/components/kernel/reference/RCBase/GeoVSurfaceShape/AnnulusDemo2.png differ
diff --git a/documentation/docs/components/kernel/reference/RCBase/GeoVSurfaceShape/Diamond.png b/documentation/docs/components/kernel/reference/RCBase/GeoVSurfaceShape/Diamond.png
new file mode 100644
index 0000000000000000000000000000000000000000..348e347132f0a3a2cf5952adb087a77ddc5c68fc
Binary files /dev/null and b/documentation/docs/components/kernel/reference/RCBase/GeoVSurfaceShape/Diamond.png differ
diff --git a/documentation/docs/components/kernel/reference/RCBase/GeoVSurfaceShape/DiamondDemo.png b/documentation/docs/components/kernel/reference/RCBase/GeoVSurfaceShape/DiamondDemo.png
new file mode 100644
index 0000000000000000000000000000000000000000..9736ca97c6e52551ab135fd35c5a912867bca645
Binary files /dev/null and b/documentation/docs/components/kernel/reference/RCBase/GeoVSurfaceShape/DiamondDemo.png differ
diff --git a/documentation/docs/components/kernel/reference/RCBase/GeoVSurfaceShape/GeoAnnulusSurface.md b/documentation/docs/components/kernel/reference/RCBase/GeoVSurfaceShape/GeoAnnulusSurface.md
new file mode 100644
index 0000000000000000000000000000000000000000..d27674bea68ec8522656c08e6324d9114e6309b6
--- /dev/null
+++ b/documentation/docs/components/kernel/reference/RCBase/GeoVSurfaceShape/GeoAnnulusSurface.md
@@ -0,0 +1,33 @@
+### GeoAnnulusSurface
+
+
+```cpp
+// === GeoAnnulusSurface === 
+
+  // Constructor:
+  // The default position of the cylinder center is at (0,0,0)
+  // Ox, Oy: the deviation of the focusing point from the default position
+  // radius_in: the inner radius of the annulus
+  // radius_out: the outer radius of the annulus
+  // phi: the span angle of the deviated circular sector, when phi = 2*PI, it is an annulus
+  GeoAnnulusSurface(double Ox, double Oy, double radius_in, double radius_out, double phi);
+
+
+  // Public Methods:
+  double getOx() const;
+  double getOy() const;
+  double getRadiusIn() const;
+  double getRadiusOut() const;
+  double getPhi() const;
+  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;
+```
+
+A `GeoAnnulusSurface` represents an annulus (or circular sector) virtual surface initialized at XOY plane. It is specified by *Ox (X-position of deviated center)*, *Oy (Y-position of deviated center)*, *Inner Radius*, *Outer Radius*, and  *Span Angle $\Phi$*, all input parameters of the constructor. The methods `getOx()`, `getOy()`, `getRadiusIn()`, `getRadiusOut()`, and `getPhi()` return the shape boundaries of the surface. When $\phi = 2\pi$, the circular section returns to an annulus.
+
+
+{{ imgutils_image_caption('RCBase/GeoVSurfaceShape/Annulus.png', 
+   alt='The GeoAnnulusSurface shape', 
+   cap='Figure: Shape boundaries of annulus surface.',
+   urlFix=False) 
+}}
diff --git a/documentation/docs/components/kernel/reference/RCBase/GeoVSurfaceShape/GeoDiamondSurface.md b/documentation/docs/components/kernel/reference/RCBase/GeoVSurfaceShape/GeoDiamondSurface.md
new file mode 100644
index 0000000000000000000000000000000000000000..7cdcc792b33673f470833d5ce7deec534bc27d8e
--- /dev/null
+++ b/documentation/docs/components/kernel/reference/RCBase/GeoVSurfaceShape/GeoDiamondSurface.md
@@ -0,0 +1,31 @@
+### GeoDiamondSurface
+
+
+```cpp
+// === GeoDiamondSurface ===
+
+  // Constructor:
+  GeoDiamondSurface (double X_bottom_half, double X_mid_half, double X_top_half, double Y_bottom_half, double Y_top_half);
+
+
+  // Public Methods:
+  double area() const;
+
+  double getXbottomHalf () const;
+  double getXmidHalf () const;
+  double getXtopHalf () const;
+  double getYbottomHalf () const;
+  double getYtopHalf () const;
+  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;
+```
+
+A `GeoDiamondSurface` represents a diamond virtual surface initialized at XOY plane. It is specified by *X half-length bottom*, *X half-length middle*, *X half-length top* in X-direction and *Y half-length bottom*, *Y half-length top* in Y-direction, which are input parameters of the constructor. The public methods `exec` and `isOnSurface` inherits methods from `GeoVSurfaceShape`, methods `getXbottomHalf()`, `getXmidHalf()`, `getXtopHalf()`, `getYbottomHalf()`, and `getYtopHalf()` return the shape boundaries of the surface. `area()` returns the area of the diamond.
+
+
+{{ imgutils_image_caption('RCBase/GeoVSurfaceShape/Diamond.png', 
+   alt='The GeoDiamondSurface shape', 
+   cap='Figure: Shape boundaries of diamond surface.',
+   urlFix=False) 
+}}
+
diff --git a/documentation/docs/components/kernel/reference/RCBase/GeoVSurfaceShape/GeoRectSurface.md b/documentation/docs/components/kernel/reference/RCBase/GeoVSurfaceShape/GeoRectSurface.md
new file mode 100644
index 0000000000000000000000000000000000000000..c225fd08347c7ca0ae02182dc45ac429ea081629
--- /dev/null
+++ b/documentation/docs/components/kernel/reference/RCBase/GeoVSurfaceShape/GeoRectSurface.md
@@ -0,0 +1,27 @@
+### GeoRectSurface
+
+
+```cpp
+// === GeoRectSurface ===
+
+   // Constructor:
+   GeoRectSurface (double XHalfLength, double YHalfLength);
+
+
+   // Public Methods:
+   double area() const;
+   double getXHalfLength () const;
+   double getYHalfLength () const;
+   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;
+```
+
+A `GeoRectSurface` represents a rectangle virtual surface initialized at the XOY plane. It is specified by a *half-length* in X-direction and a *half-length* in Y-direction, which are input parameters of the constructor.  Methods `getXHalfLength()`, `getYHalfLength()` return the shape boundaries of the surface. `area()` returns the area.
+
+
+{{ imgutils_image_caption('RCBase/GeoVSurfaceShape/Rect.png', 
+   alt='The GeoRectSurface shape', 
+   cap='Figure: Shape boundaries of rectangle surface.',
+   urlFix=False) 
+}}
+
diff --git a/documentation/docs/components/kernel/reference/RCBase/GeoVSurfaceShape/GeoTrapezoidSurface.md b/documentation/docs/components/kernel/reference/RCBase/GeoVSurfaceShape/GeoTrapezoidSurface.md
new file mode 100644
index 0000000000000000000000000000000000000000..0f8793bd50fb3247bc98278530976b732e2dad7b
--- /dev/null
+++ b/documentation/docs/components/kernel/reference/RCBase/GeoVSurfaceShape/GeoTrapezoidSurface.md
@@ -0,0 +1,27 @@
+### GeoTrapezoidSurface
+
+
+```cpp
+// === GeoTrapezoidSurface ===
+
+  // Constructor:
+  GeoTrapezoidSurface (double XHalfLengthMin, double XHalfLengthMax, double YHalfLength);
+
+
+  // Public Methods:
+  double area() const;
+  double getXHalfLengthMin () const;
+  double getXHalfLengthMax () const;
+  double getYHalfLength () const;
+  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;
+```
+
+A `GeoTrapezoidSurface` represents a trapezoid virtual surface initialized at XOY plane. It is specified by *X half-length minimum* (or bottom half-length) and *X half-length maximum* (or top half-length) in X-direction and *Y half-length* in Y-direction, which are input parameters of the constructor. The methods `getXHalfLengthMin()`, `getXHalfLengthMax`, and `getYHalfLength()` return the shape boundaries of the surface. `area()` returns the area of the trapezoid.
+
+
+{{ imgutils_image_caption('RCBase/GeoVSurfaceShape/Trapezoid.png', 
+   alt='The GeoTrapezoidSurface shape', 
+   cap='Figure: Shape boundaries of trapezoid surface.',
+   urlFix=False) 
+}}
diff --git a/documentation/docs/components/kernel/reference/RCBase/GeoVSurfaceShape/GeoVSurface.md b/documentation/docs/components/kernel/reference/RCBase/GeoVSurfaceShape/GeoVSurface.md
new file mode 100644
index 0000000000000000000000000000000000000000..fc766ba8e2bae75b5783b08af9ccfe4e4338db84
--- /dev/null
+++ b/documentation/docs/components/kernel/reference/RCBase/GeoVSurfaceShape/GeoVSurface.md
@@ -0,0 +1,23 @@
+### GeoVSurface
+
+
+The class GeoVSurface takes [Virtual Surface Shape](#geovsurfaceshape) as input to the constructor, and is used to embed the shape within the geometry tree.  It provides access to the SurfaceShape and can determine whether a triplet of coordinates lies on and within the surface. 
+
+
+```cpp
+// === GeoVSurface ===
+
+   //Constructor:
+   GeoVSurface(const GeoVSurfaceShape *SurfShape);
+
+   //Public Methods:
+   const GeoVSurfaceShape* getShape () const
+   virtual void exec(GeoNodeAction *action) const;
+   bool isOnSurface(const double Px, const double Py, const double Pz) const;
+```
+For example, if we want to construct a rectangle virtual surface node, we just need to first construct the surface shape and then construct the virtual surface node:
+
+```cpp
+   GeoRectSurface* rectangle = new GeoRectSurface(5, 9); // 5 is X half length and 9 is Y half length
+   GeoVSurface* surface = new GeoVSurface(rectangle);    // construct a virtual surface node in rectangle shape   
+```
\ No newline at end of file
diff --git a/documentation/docs/components/kernel/reference/RCBase/GeoVSurfaceShape/GeoVSurfaceShape.md b/documentation/docs/components/kernel/reference/RCBase/GeoVSurfaceShape/GeoVSurfaceShape.md
new file mode 100644
index 0000000000000000000000000000000000000000..46ed45d46d2574464124522c2ed690b42d27b151
--- /dev/null
+++ b/documentation/docs/components/kernel/reference/RCBase/GeoVSurfaceShape/GeoVSurfaceShape.md
@@ -0,0 +1,19 @@
+### GeoVSurfaceShape
+
+
+The abstract class GeoVSurfaceShape defines the interface for subclasses, specifically the `isOnSurface` method. 
+
+```cpp
+// === GeoVSurfaceShape ===
+
+   //Constructor:
+   GeoVSurfaceShape () = default;
+
+   //Public Methods:
+   //Executes a GeoShapeAction on the surface shape
+   virtual void exec (GeoShapeAction *action) const = 0;
+   //Determines whether the point(Px,Py,Pz) is on the surface, can be overriden by child methods
+   virtual bool isOnSurface (const double Px, const double Py, const double Pz, const GeoTrf::Transform3D & trans) const = 0;
+```
+
+All surface shapes are initialized at X-O-Y plane. In *gmex* when visualized, the side initially facing at +Z direction is in red color, while the side initially facing at -Z direction is in green color. The following subsections show specific examples of virtual surface shapes.
\ No newline at end of file
diff --git a/documentation/docs/components/kernel/reference/RCBase/GeoVSurfaceShape/Rect.png b/documentation/docs/components/kernel/reference/RCBase/GeoVSurfaceShape/Rect.png
new file mode 100644
index 0000000000000000000000000000000000000000..21ae1c57a0d18f6df62f9ba97ab1da915052c3c6
Binary files /dev/null and b/documentation/docs/components/kernel/reference/RCBase/GeoVSurfaceShape/Rect.png differ
diff --git a/documentation/docs/components/kernel/reference/RCBase/GeoVSurfaceShape/RectDemo.png b/documentation/docs/components/kernel/reference/RCBase/GeoVSurfaceShape/RectDemo.png
new file mode 100644
index 0000000000000000000000000000000000000000..82796f1e0c82c7bfe8887df86da60fbf08892139
Binary files /dev/null and b/documentation/docs/components/kernel/reference/RCBase/GeoVSurfaceShape/RectDemo.png differ
diff --git a/documentation/docs/components/kernel/reference/RCBase/GeoVSurfaceShape/SurfaceClassRef.png b/documentation/docs/components/kernel/reference/RCBase/GeoVSurfaceShape/SurfaceClassRef.png
new file mode 100644
index 0000000000000000000000000000000000000000..9d42ba7aa6a92752c85ca076975e0904e2548175
Binary files /dev/null and b/documentation/docs/components/kernel/reference/RCBase/GeoVSurfaceShape/SurfaceClassRef.png differ
diff --git a/documentation/docs/components/kernel/reference/RCBase/GeoVSurfaceShape/Trapezoid.png b/documentation/docs/components/kernel/reference/RCBase/GeoVSurfaceShape/Trapezoid.png
new file mode 100644
index 0000000000000000000000000000000000000000..deb177c9ecf1fc22e06efa7fe06f338ac7b0d25b
Binary files /dev/null and b/documentation/docs/components/kernel/reference/RCBase/GeoVSurfaceShape/Trapezoid.png differ
diff --git a/documentation/docs/components/kernel/reference/RCBase/GeoVSurfaceShape/TrapezoidDemo.png b/documentation/docs/components/kernel/reference/RCBase/GeoVSurfaceShape/TrapezoidDemo.png
new file mode 100644
index 0000000000000000000000000000000000000000..da24d089cece86b0b68ddb735bbea22238fddd7c
Binary files /dev/null and b/documentation/docs/components/kernel/reference/RCBase/GeoVSurfaceShape/TrapezoidDemo.png differ
diff --git a/documentation/docs/components/kernel/reference/VirtualSurface.md b/documentation/docs/components/kernel/reference/VirtualSurface.md
new file mode 100644
index 0000000000000000000000000000000000000000..d2c2b184e7cd52616973aa1ea7e8ad19a9ab753e
--- /dev/null
+++ b/documentation/docs/components/kernel/reference/VirtualSurface.md
@@ -0,0 +1,14 @@
+
+## Virtual Surfaces
+
+### Introduction
+
+Virtual surface classes are light weight classes describing geometrical objects devoid of any material properties. They are introduced for particle tracking and fast simulation purposes. Like `GeoVFullPhysVol`, the `GeoVSurface` class caches its absolute position locally and responds to alignment updates further up the tree. Virtual surface nodes may be parented within *physical volumes*, but they do not have any child nodes of any kind. The node actions are adapted to work with surface classes as well.
+
+{{ imgutils_image_caption('RCBase/GeoVSurfaceShape/SurfaceClassRef.png', 
+   alt='The GeoVSurface inheritance', 
+   cap='Figure: GeoVSurface inheritance relations.',
+   urlFix=False) 
+}}
+
+{% include 'components/kernel/reference/RCBase/GeoVSurfaceShape/GeoVSurface.md' %}
\ No newline at end of file
diff --git a/documentation/docs/components/kernel/reference/VirtualSurfaceShape.md b/documentation/docs/components/kernel/reference/VirtualSurfaceShape.md
new file mode 100644
index 0000000000000000000000000000000000000000..3eca9cdd3192f1a9bc162a2025241612e10f86dd
--- /dev/null
+++ b/documentation/docs/components/kernel/reference/VirtualSurfaceShape.md
@@ -0,0 +1,35 @@
+
+
+##	Virtual Surface Shape
+
+
+### Introduction
+
+The virtual surface shape classes are designed to describe the geometry of two dimensional virtual surfaces. The specific surface shape classes are used as input variables of `GeoVSurface`, as described in section [GeoVSurface](#geovsurface). The set of surfaces is extensible.
+
+
+
+In the table here below, a list of the existing virtual surface shapes in the GeoModelKernel package is given:
+
+
+| Class   | Virtual Surface Shape |
+| ------- | ----- |
+| [GeoRectSurface](#georectsurface)  | Rectangle |
+| [GeoTrapezoidSurface](#geotrapezoidsurface) | Trapezoid |
+| [GeoDiamondSurface](#geodiamondsurface) | Diamond |
+| [GeoAnnulusSurface](#geoannulussurface) | Annulus / Circular Sector |
+
+The surface shapes can be visualized by *gmex*, and some screenshots appear in the following documentation.
+
+
+{% include 'components/kernel/reference/RCBase/GeoVSurfaceShape/GeoVSurfaceShape.md' %}
+
+{% include 'components/kernel/reference/RCBase/GeoVSurfaceShape/GeoRectSurface.md' %}
+
+{% include 'components/kernel/reference/RCBase/GeoVSurfaceShape/GeoTrapezoidSurface.md' %}
+
+{% include 'components/kernel/reference/RCBase/GeoVSurfaceShape/GeoDiamondSurface.md' %}
+
+{% include 'components/kernel/reference/RCBase/GeoVSurfaceShape/GeoAnnulusSurface.md' %}
+
+
diff --git a/documentation/docs/components/kernel/reference/index.md b/documentation/docs/components/kernel/reference/index.md
index 3914f0ebd38ec659685ce16fce1d9b856323b4c1..2289376fa8fa2231c0d4ecdf6470dcaedf8bcded 100644
--- a/documentation/docs/components/kernel/reference/index.md
+++ b/documentation/docs/components/kernel/reference/index.md
@@ -15,6 +15,10 @@ Detailed descriptions of the geometry kernel classes follow.
 
 {% include 'components/kernel/reference/LogicalVolumes.md' %}
 
+{% include 'components/kernel/reference/VirtualSurfaceShape.md' %}
+
+{% include 'components/kernel/reference/VirtualSurface.md' %}
+
 {% include 'components/kernel/reference/PhysicalVolumes.md' %}
 
 {% include 'components/kernel/reference/Transformations.md' %}