diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoBox.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoBox.h
index b024be9ad71ee38aade79b6114fda765cc655bfb..d4ed518df9493bc04ec3251a92dbb70fb7bb1add 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoBox.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoBox.h
@@ -10,34 +10,41 @@
 class GeoBox : public GeoShape
 {
  public:
-  //	Constructor for the BOX.
+  //    Constructor for the BOX
   GeoBox (double XHalfLength, double YHalfLength, double ZHalfLength);
 
-  //	Returns the volume of the shape, for mass inventory
+  //    Returns the volume of the shape, for mass inventory
   virtual double volume () const;
 
-  //	Returns the BOX shape type, as a string.
+  //    Returns the bonding box of the shape
+  virtual void extent (double& xmin, double& ymin, double& zmin,
+                       double& xmax, double& ymax, double& zmax) const;
+
+  //    Returns true if the shape contains the point, false otherwise.
+  virtual bool contains (double x, double y, double z) const;
+
+  //    Returns the BOX shape type, as a string
   virtual const std::string & type () const;
 
-  //	Returns the BOX shape type, as a coded integer.
+  //    Returns the BOX shape type, as a coded integer
   virtual ShapeType typeID () const;
 
-  //	 Executes a GeoShapeAction.
+  //     Executes a GeoShapeAction
   virtual void exec (GeoShapeAction *action) const;
 
-  //	For type identification.
+  //    For type identification
   static const std::string& getClassType ();
 
-  //	For type identification.,
+  //    For type identification
   static ShapeType getClassTypeID ();
 
-  //	Half length in the x-direction.
+  //    Half length in the x-direction
   const double& getXHalfLength () const;
 
-  //	Half-length in the y direction.
+  //    Half-length in the y direction
   const double& getYHalfLength () const;
 
-  //	Half-length in the z direction.
+  //    Half-length in the z direction
   const double& getZHalfLength () const;
 
  protected:
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoCons.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoCons.h
index 052e59d381bf07a077d7b70dc0580a4aeed28510..58528bb868f278f58ac46651f359623b98754696 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoCons.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoCons.h
@@ -12,52 +12,59 @@ class GeoCons : public GeoShape
  public:
   GeoCons (double RMin1, double RMin2, double RMax1, double RMax2, double DZ, double SPhi, double DPhi);
 
-  //	Returns the volume of the shape, for mass inventory
+  //    Returns the volume of the shape, for mass inventory.
   virtual double volume () const;
 
-  //	Returns the CONS shape type, as a string.
+  //    Returns the bonding box of the shape.
+  virtual void extent (double& xmin, double& ymin, double& zmin,
+                       double& xmax, double& ymax, double& zmax) const;
+
+  //    Returns true if the shape contains the point, false otherwise.
+  virtual bool contains (double x, double y, double z) const;
+
+  //    Returns the CONS shape type, as a string.
   virtual const std::string & type () const;
-  
-  //	Returns the CONS shape type, as a coded integer.
+
+  //    Returns the CONS shape type, as a coded integer.
   virtual ShapeType typeID () const;
-  
-  //	Executes a GeoShapeAction
+
+  //    Executes a GeoShapeAction.
   virtual void exec (GeoShapeAction *action) const;
-  
-  //	For type identification.
+
+  //    For type identification.
   static const std::string& getClassType ();
-  
-  //	For type identification.
+
+  //    For type identification.
   static ShapeType getClassTypeID ();
-  
-  //	Inside radius at -dZ
+
+  //    Inside radius at -dZ
   const double& getRMin1 () const;
-  
-  //	Inside radius at +dZ
+
+  //    Inside radius at +dZ
   const double& getRMin2 () const;
-  
-  //	Outside radius at -dZ
+
+  //    Outside radius at -dZ
   const double& getRMax1 () const;
-  
-  //	Outside radius at +dZ
+
+  //    Outside radius at +dZ
   const double& getRMax2 () const;
-  
-  //	Half length in Z direction.
+
+  //    Half length in Z direction.
   const double& getDZ () const;
-  
-  //	Starting angle of the segment in radians
+
+  //    Starting angle of the segment in radians.
   const double& getSPhi () const;
-  
-  //	Delta angle of the segment in radians.
+
+  //    Delta angle of the segment in radians.
   const double& getDPhi () const;
-  
+
  protected:
   virtual ~GeoCons();
-  
+
  private:
   GeoCons(const GeoCons &right);
   GeoCons & operator=(const GeoCons &right);
-  
+
   static const std::string s_classType;
   static const ShapeType s_classTypeID;
 
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoEllipticalTube.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoEllipticalTube.h
index a374b8e6db7372c8d8c063810c27cbdc97b8acaa..953c8405e05ba42c4596bc5f8c5f008da571bab1 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoEllipticalTube.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoEllipticalTube.h
@@ -7,7 +7,7 @@
 
 /**
  * @class: GeoEllipticalTube
- *  
+ *
  * @brief This is a tube with elliptical cross section
  * The equation of the surface in x/y is 1.0 = (x/dx)**2 + (y/dy)**2
  */
@@ -17,27 +17,48 @@
 class GeoEllipticalTube : public GeoShape
 {
  public:
-  GeoEllipticalTube(double XHalfLength, double YHalfLength, double ZHalfLength); 
+  //    Constructor for the ELLIPTICAL TUBE
+  GeoEllipticalTube(double XHalfLength, double YHalfLength, double ZHalfLength);
 
+  //    Returns the volume of the shape, for mass inventory
   virtual double volume () const;
 
+  //    Returns the bonding box of the shape
+  virtual void extent (double& xmin, double& ymin, double& zmin,
+                       double& xmax, double& ymax, double& zmax) const;
+
+  //    Returns true if the shape contains the point, false otherwise.
+  virtual bool contains (double x, double y, double z) const;
+
+  //    Returns the ELLIPTICAL TUBE shape type, as a string
   virtual const std::string & type () const;
+
+  //    Returns the ELLIPTICAL TUBE 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 ();
 
+  //    X semi-axis
   const double& getXHalfLength() const;
+
+  //    Y semi-axis
   const double& getYHalfLength() const;
+
+  //    Elliptical tube half-length in the Z direction
   const double& getZHalfLength() const;
 
  protected:
   virtual ~GeoEllipticalTube();
 
  private:
-  
+
   GeoEllipticalTube(const GeoEllipticalTube &right);
   GeoEllipticalTube & operator=(const GeoEllipticalTube &right);
 
@@ -49,7 +70,6 @@ class GeoEllipticalTube : public GeoShape
   double m_zHalfLength;
 };
 
-
 inline const std::string& GeoEllipticalTube::getClassType ()
 {
   return s_classType;
@@ -75,5 +95,4 @@ inline const double& GeoEllipticalTube::getZHalfLength() const
   return m_zHalfLength;
 }
 
-
 #endif
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoGenericTrap.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoGenericTrap.h
index 7db180074407007f0d2951ab6be2cb6288836dca..b8a310ad201e6abf9bd97346ad32eadea045d622 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoGenericTrap.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoGenericTrap.h
@@ -17,17 +17,35 @@ class GeoGenericTrap : public GeoShape
  public:
   GeoGenericTrap(double ZHalfLength, const GeoGenericTrapVertices& Vertices);
 
+  //    Returns the volume of the shape, for mass inventory.
   virtual double volume() const;
 
+  //    Returns the bonding box of the shape.
+  virtual void extent (double& xmin, double& ymin, double& zmin,
+                       double& xmax, double& ymax, double& zmax) const;
+
+  //    Returns true if the shape contains the point, false otherwise.
+  virtual bool contains (double x, double y, double z) const;
+
+  //    Returns the GENERIC TRAP shape type, as a string.
   virtual const std::string& type() const;
+
+  //    Returns the GENERIC TRAP shape type, as a coded integer.
   virtual ShapeType typeID() const;
 
+  //    For type identification.
   static const std::string& getClassType();
+
+  //    For type identification.
   static ShapeType getClassTypeID();
 
+  //    Executes a GeoShapeAction.
   virtual void exec(GeoShapeAction *action) const;
 
+  //    Z half length.
   double getZHalfLength() const;
+
+  //    Vector of vertices.
   const GeoGenericTrapVertices& getVertices() const;
 
  protected:
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoPara.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoPara.h
index 7072aee43ce53ba6fd65af0aac16541152363b96..3b7485b6856763e33b69b01e82b585d89e8ab072 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoPara.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoPara.h
@@ -10,53 +10,60 @@
 class GeoPara : public GeoShape
 {
  public:
-  //	Constructor for the BOX.
+  //    Constructor for the BOX.
   GeoPara (double XHalfLength, double YHalfLength, double ZHalfLength, double Alpha, double Theta, double Phi);
-  
-  //	Returns the volume of the shape, for mass inventory
+
+  //    Returns the volume of the shape, for mass inventory.
   virtual double volume () const;
-  
-  //	Returns the PARA shape type, as a string.
+
+  //    Returns the bonding box of the shape.
+  virtual void extent (double& xmin, double& ymin, double& zmin,
+                       double& xmax, double& ymax, double& zmax) const;
+
+  //    Returns true if the shape contains the point, false otherwise.
+  virtual bool contains (double x, double y, double z) const;
+
+  //    Returns the PARA shape type, as a string.
   virtual const std::string & type () const;
-  
-  //	Returns the PARA shape type, as a coded integer.
+
+  //    Returns the PARA shape type, as a coded integer.
   virtual ShapeType typeID () const;
-  
-  //	Executes a GeoShapeAction
+
+  //    Executes a GeoShapeAction
   virtual void exec (GeoShapeAction *action) const;
-  
-  //	For type identification.
+
+  //    For type identification.
   static const std::string& getClassType ();
-  
-  //	For type identification.
+
+  //    For type identification.
   static ShapeType getClassTypeID ();
-  
-  //	Half length in the x-direction.
+
+  //    Half length in the x-direction.
   const double& getXHalfLength () const;
-  
-  //	Half-length in the y direction.
+
+  //    Half-length in the y direction.
   const double& getYHalfLength () const;
-  
-  //	Half-length in the z direction.
+
+  //    Half-length in the z direction.
   const double& getZHalfLength () const;
-  
-  //	Polar (theta) angle.
+
+  //    Polar (theta) angle.
   const double& getTheta () const;
-  
-  //	The angle alpha...between the two sides of the top face
-  //	of the parallelapiped.
+
+  //    The angle alpha...between the two sides of the top face
+  //    of the parallelapiped.
   const double& getAlpha () const;
-  
-  //	Azimuthal (phi) angle.
+
+  //    Azimuthal (phi) angle.
   const double& getPhi () const;
-  
+
  protected:
   virtual ~GeoPara();
-  
+
  private:
   GeoPara(const GeoPara &right);
   GeoPara & operator=(const GeoPara &right);
-  
+
   static const std::string s_classType;
   static const ShapeType s_classTypeID;
 
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoPcon.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoPcon.h
index c08797c8d8fefca975fd148c8d2089663adbe21a..90979919813c513e96d03166c29dcf9fc5cce5d2 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoPcon.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoPcon.h
@@ -10,7 +10,7 @@
  *
  * @brief This shape represents a polycone.
  * Specify the starting angle and delta angle (subtended angle)
- * of the polycone first, Then add at least two planes with the
+ * of the polycone first, then add at least two planes with the
  * addPlane( double zPlane, double rInner,  double rOuter) method.
  */
 
@@ -20,76 +20,83 @@
 class GeoPcon : public GeoShape
 {
  public:
-  //	Constructor for the PCON.  Note that the constructor
-  //	does not fully build this object.  The PCON is not valid
-  //	until at least two polygon planes have been added.
+  //    Constructor for the PCON.  Note that the constructor
+  //    does not fully build this object.  The PCON is not valid
+  //    until at least two polygon planes have been added.
   GeoPcon (double SPhi, double DPhi);
 
-  //	Returns the volume of the shape, for mass inventory
+  //    Returns the volume of the shape, for mass inventory
   virtual double volume () const;
-  
-  //	Returns the PCON shape type, as a string.
+
+  //    Returns the bonding box of the shape
+  virtual void extent (double& xmin, double& ymin, double& zmin,
+                       double& xmax, double& ymax, double& zmax) const;
+
+  //    Returns true if the shape contains the point, false otherwise.
+  virtual bool contains (double x, double y, double z) const;
+
+  //    Returns the PCON shape type, as a string.
   virtual const std::string & type () const;
-  
-  //	Returns the PCON shape type, as a coded integer.
+
+  //    Returns the PCON shape type, as a coded integer.
   virtual ShapeType typeID () const;
-  
-  //	Add another plane to the polycone  A minimum of two
-  //	planes are required to create a valid polycone.
+
+  //    Add another plane to the polycone  A minimum of two
+  //    planes are required to create a valid polycone.
   void addPlane (double ZPlane, double RMinPlane, double RMaxPlane);
-  
-  //	Returns the number of planes that have been created for
-  //	the polycone.
+
+  //    Returns the number of planes that have been created for
+  //    the polycone.
   unsigned int getNPlanes () const;
-  
-  //	True if the polycone has at least two planes.  False
-  //	otherwise.
+
+  //    True if the polycone has at least two planes. False
+  //    otherwise.
   bool isValid () const;
-  
-  //	Get the Z Position of the specified plane.
+
+  //    Get the Z Position of the specified plane.
   const double & getZPlane (unsigned int i) const;
-  
-  //	Get the RMin of the specified plane.
+
+  //    Get the RMin of the specified plane.
   const double & getRMinPlane (unsigned int i) const;
-  
-  //	Get the Z Position of the specified plane.
+
+  //    Get the Z Position of the specified plane.
   const double & getRMaxPlane (unsigned int i) const;
-  
-  //	Executes a GeoShapeAction
+
+  //    Executes a GeoShapeAction
   virtual void exec (GeoShapeAction *action) const;
-  
-  //	For type identification.
+
+  //    For type identification.
   static const std::string& getClassType ();
-  
-  //	For type identification.
+
+  //    For type identification.
   static ShapeType getClassTypeID ();
-  
-  //	Starting angle of the segment in radians.
+
+  //    Starting angle of the segment in radians.
   const double& getSPhi () const;
-  
-  //	Delta angle of the segment in radians.
+
+  //    Delta angle of the segment in radians.
   const double& getDPhi () const;
 
  protected:
   virtual ~GeoPcon();
-  
+
  private:
   GeoPcon(const GeoPcon &right);
   GeoPcon & operator=(const GeoPcon &right);
-  
+
   static const std::string s_classType;
   static const ShapeType s_classTypeID;
 
   double m_sPhi;
   double m_dPhi;
 
-  //	Z Position of poly-cone planes.
+  //    Z Position of poly-cone planes.
   std::vector<double> m_zPlane;
-  
-  //	Minimum radius of poly-cone planes.
+
+  //    Minimum radius of poly-cone planes.
   std::vector<double> m_rMinPlane;
 
-  //	Maximum radius of poly-cone planes.
+  //    Maximum radius of poly-cone planes.
   std::vector<double> m_rMaxPlane;
 };
 
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoPgon.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoPgon.h
index 12f22cb28afea3254618d157e9a12df7d5b129ca..ef179c32701d7df43806216d855a66d68417f9d9 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoPgon.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoPgon.h
@@ -11,66 +11,73 @@
 class GeoPgon : public GeoShape
 {
  public:
-  //	Constructor for the PGON.  Note that the constructor
-  //	does not fully build this object.  The PGON is not valid
-  //	until at least two polygon planes have been added.
+  //    Constructor for the PGON.  Note that the constructor
+  //    does not fully build this object.  The PGON is not valid
+  //    until at least two polygon planes have been added.
   GeoPgon (double SPhi, double DPhi, unsigned int NSides);
 
-  //	Returns the volume of the shape, for mass inventory
+  //    Returns the volume of the shape, for mass inventory
   virtual double volume () const;
-  
-  //	Returns the PGON shape type, as a string.
+
+  //    Returns the bonding box of the shape
+  virtual void extent (double& xmin, double& ymin, double& zmin,
+                       double& xmax, double& ymax, double& zmax) const;
+
+  //    Returns true if the shape contains the point, false otherwise.
+  virtual bool contains (double x, double y, double z) const;
+
+  //    Returns the PGON shape type, as a string.
   virtual const std::string & type () const;
-  
-  //	Returns the PGON shape type, as a coded integer.
+
+  //    Returns the PGON shape type, as a coded integer.
   virtual ShapeType typeID () const;
-  
-  //	Add another plane to the polygon.  A minimum of two
-  //	planes are required to create a valid polygon.
+
+  //    Add another plane to the polygon. A minimum of two
+  //    planes are required to create a valid polygon.
   void addPlane (double ZPlane, double RMinPlane, double RMaxPlane);
-  
-  //	Returns the number of planes that have been created for
-  //	the polygon.
+
+  //    Returns the number of planes that have been created for
+  //    the polygon.
   unsigned int getNPlanes () const;
-  
-  //	True if the polygon has at least two planes.  False
-  //	otherwise.
+
+  //    True if the polygon has at least two planes.  False
+  //    otherwise.
   bool isValid () const;
-  
-  //	Get the Z Position of the specified plane.
+
+  //    Get the Z Position of the specified plane.
   const double & getZPlane (unsigned int i) const;
-  
-  //	Get the RMin of the specified plane.
+
+  //    Get the RMin of the specified plane.
   const double & getRMinPlane (unsigned int i) const;
-  
-  //	Get the Z Position of the specified plane.
+
+  //    Get the Z Position of the specified plane.
   const double & getRMaxPlane (unsigned int i) const;
-  
-  //	Executes a GeoShapeAction
+
+  //    Executes a GeoShapeAction
   virtual void exec (GeoShapeAction *action) const;
-  
-  //	For type identification.
+
+  //    For type identification.
   static const std::string& getClassType ();
-  
-  //	For type identification.
+
+  //    For type identification.
   static ShapeType getClassTypeID ();
-  
-  //	Starting angle of the segment in radians.
+
+  //    Starting angle of the segment in radians.
   const double& getSPhi () const;
-  
-  //	Delta angle of the segment in radians.
+
+  //    Delta angle of the segment in radians.
   const double& getDPhi () const;
-  
-  //	Number of sides in each polygonal segment.
+
+  //    Number of sides in each polygonal segment.
   const unsigned int& getNSides () const;
-  
+
  protected:
   virtual ~GeoPgon();
-  
+
  private:
   GeoPgon(const GeoPgon &right);
   GeoPgon & operator=(const GeoPgon &right);
-  
+
   static const std::string s_classType;
   static const ShapeType s_classTypeID;
 
@@ -78,13 +85,13 @@ class GeoPgon : public GeoShape
   double m_dPhi;
   unsigned int m_nSides;
 
-  //	Z Position of polygon planes.
+  //    Z Position of polygon planes.
   std::vector<double> m_zPlane;
 
-  //	Minimum radius of polygon planes.
+  //    Minimum radius of polygon planes.
   std::vector<double> m_rMinPlane;
 
-  //	Maximum radius of polygon planes.
+  //    Maximum radius of polygon planes.
   std::vector<double> m_rMaxPlane;
 };
 
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoShape.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoShape.h
index 8241e30b78610c63b8d066296dfdb009f0bfe93c..fc3b4387c65e8b7f96ea0bd3a06ba912374665d6 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoShape.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoShape.h
@@ -16,12 +16,12 @@
  *           if (myShape->typeId()==GeoBox::classTypeId()) {
  *             .....
  *           }
- *	This avoids the need for dynamic casting with these
- *	classes.
+ *      This avoids the need for dynamic casting with these
+ *      classes.
  *
- *	Also noteworthy: Shapes are allocated on the heap and
- *	deleted automatically when their reference count falls
- *	to zero.
+ *      Also noteworthy: Shapes are allocated on the heap and
+ *      deleted automatically when their reference count falls
+ *      to zero.
  */
 
 #include "GeoModelKernel/RCBase.h"
@@ -38,36 +38,48 @@ class GeoShapeAction;
 class GeoShape : public RCBase
 {
  public:
-  // Constructor for shape.  Must provide the name, a string to identify this shape.
+  //    Constructor for shape. Must provide the name, a string to identify this shape.
   GeoShape ();
 
-  //	Returns the volume of the shape, for mass inventory
-  virtual double volume () const = 0;
+  //    Returns the volume of the shape, for mass inventory.
+  virtual double volume () const;
 
-  //	Boolean OR operation for shapes
+  //    Returns the bonding box of the shape.
+  virtual void extent (double& xmin, double& ymin, double& zmin,
+                       double& xmax, double& ymax, double& zmax) const = 0;
+
+  //    Returns true if the shape contains the point, false otherwise.
+  virtual bool contains (double x, double y, double z) const = 0;
+
+  //    Boolean OR operation for shapes.
   const GeoShapeUnion & add (const GeoShape& shape) const;
-  
-  //	Subtraction operation for shapes.
+
+  //    Subtraction operation for shapes.
   const GeoShapeSubtraction & subtract (const GeoShape& shape) const;
-  
-  //	Intersection of shapes.
+
+  //    Intersection of shapes.
   const GeoShapeIntersection & intersect (const GeoShape& shape) const;
-  
-  //	Shift shapes, for boolean operations.
+
+  //    Shift shapes, for boolean operations.
   const GeoShapeShift & operator << (const GeoTrf::Transform3D &shift) const;
-  
-  //	Returns the shape type, as a string.
+
+  //    Returns the shape type, as a string.
   virtual const std::string & type () const = 0;
-  
-  //	Returns the shape type, as an coded integer.
+
+  //    Returns the shape type, as an coded integer.
   virtual ShapeType typeID () const = 0;
-  
-  //	Executes a GeoShapeAction
+
+  //    Executes a GeoShapeAction
   virtual void exec (GeoShapeAction *action) const = 0;
-  
+
  protected:
   virtual ~GeoShape();
 
+  //    Returns the bounding box of the specified disk. This method is used
+  //    for calculation of the extend of a tube, cone, polycone, torus.
+  static void diskExtent(double rmin, double rmax, double sphi, double dphi,
+                         double& xmin, double& ymin, double& xmax, double& ymax);
+
  private:
   GeoShape(const GeoShape &right);
   GeoShape & operator=(const GeoShape &right);
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoShapeIntersection.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoShapeIntersection.h
index f01bb3dcec726b0a70bca5409134c8f7dab68173..6d9ee86eb22f684074cd8b762405c18454867785 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoShapeIntersection.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoShapeIntersection.h
@@ -7,43 +7,45 @@
 
 #include "GeoModelKernel/GeoShape.h"
 
-
 #ifndef _GeoShapePersistification_On_
   class Persistifier;
 #endif
 
-
 class GeoShapeIntersection : public GeoShape
 {
  public:
-
-
-
-  //	Constructor taking two shape operands.
+  //    Constructor taking two shape operands
   GeoShapeIntersection (const GeoShape* A, const GeoShape* B);
 
-  //	Returns the volume of the shape, for mass inventory
+  //    Returns the volume of the shape, for mass inventory
   virtual double volume () const;
 
-  //	Returns the AND shape type, as a string.
+  //    Returns the bonding box of the shape
+  virtual void extent (double& xmin, double& ymin, double& zmin,
+                       double& xmax, double& ymax, double& zmax) const;
+
+  //    Returns true if the shape contains the point, false otherwise
+  virtual bool contains (double x, double y, double z) const;
+
+  //    Returns the AND shape type, as a string.
   virtual const std::string & type () const;
 
-  //	Returns the AND shape type, as a coded integer.
+  //    Returns the AND shape type, as a coded integer.
   virtual ShapeType typeID () const;
 
-  //	Returns the first operand being ANDed
+  //    Returns the first operand being ANDed
   const GeoShape* getOpA () const;
 
-  //	Returns the second operand being ANDed.
+  //    Returns the second operand being ANDed.
   const GeoShape* getOpB () const;
 
-  //	Executes a GeoShapeAction
+  //    Executes a GeoShapeAction
   virtual void exec (GeoShapeAction *action) const;
 
-  //	For type identification.
+  //    For type identification.
   static const std::string& getClassType ();
 
-  //	For type identification.
+  //    For type identification.
   static ShapeType getClassTypeID ();
 
  protected:
@@ -53,12 +55,15 @@ class GeoShapeIntersection : public GeoShape
   GeoShapeIntersection(const GeoShapeIntersection &right);
   GeoShapeIntersection & operator=(const GeoShapeIntersection &right);
 
-  //	The first shape operand in the AND operation.
+  //    The first shape operand in the AND operation.
   const GeoShape* m_opA;
 
-  //	The second shape operand in the AND operation.
+  //    The second shape operand in the AND operation.
   const GeoShape* m_opB;
 
+  //    Cached volume
+  mutable double fVolume = -1.;
+
   static const std::string s_classType;
   static const ShapeType s_classTypeID;
 
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoShapeShift.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoShapeShift.h
index 870bf304d9ec48be8e690ab8de7cb4e2338bde4b..957b6f9ef847f653aae07ef566c8e18d3c3d73a8 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoShapeShift.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoShapeShift.h
@@ -17,30 +17,38 @@
 class GeoShapeShift : public GeoShape
 {
  public:
+  //    Constructor
   GeoShapeShift (const GeoShape* A, const GeoTrf::Transform3D &X);
 
-  //	Returns the volume of the shape, for mass inventory
+  //    Returns the volume of the shape, for mass inventory
   virtual double volume () const;
 
-  //	Returns the OR shape type, as a string.
+  //    Returns the bonding box of the shape
+  virtual void extent (double& xmin, double& ymin, double& zmin,
+                       double& xmax, double& ymax, double& zmax) const;
+
+  //    Returns true if the shape contains the point, false otherwise
+  virtual bool contains (double x, double y, double z) const;
+
+  //    Returns the OR shape type, as a string.
   virtual const std::string & type () const;
 
-  //	Returns the OR shape type, as a coded integer.
+  //    Returns the OR shape type, as a coded integer.
   virtual ShapeType typeID () const;
 
-  //	Returns the first operand being ORed
+  //    Returns the first operand being ORed
   const GeoShape* getOp () const;
 
-  //	Returns the shift of this shape.
+  //    Returns the shift of this shape.
   const GeoTrf::Transform3D & getX () const;
 
-  //	Executes a GeoShapeAction
+  //    Executes a GeoShapeAction
   virtual void exec (GeoShapeAction *action) const;
 
-  //	For type identification.
+  //    For type identification.
   static const std::string& getClassType ();
 
-  //	For type identification.
+  //    For type identification.
   static ShapeType getClassTypeID ();
 
  protected:
@@ -50,10 +58,10 @@ class GeoShapeShift : public GeoShape
   GeoShapeShift(const GeoShapeShift &right);
   GeoShapeShift & operator=(const GeoShapeShift &right);
 
-  //	The shape operand in the NOT operation.
+  //    The shape operand in the NOT operation.
   const GeoShape* m_op;
 
-  //	Gives the amount by which the volume is shifted.
+  //    Gives the amount by which the volume is shifted.
   GeoTrf::Transform3D m_shift;
 
   static const std::string s_classType;
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoShapeSubtraction.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoShapeSubtraction.h
index 864112011a1f222c3a186d9f8ab34435875b75d7..66fb34aca5890b0870f8f7af834c6d69ead435e0 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoShapeSubtraction.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoShapeSubtraction.h
@@ -16,30 +16,38 @@
 class GeoShapeSubtraction : public GeoShape
 {
  public:
+  //    Constructor taking two shape operands
   GeoShapeSubtraction (const GeoShape* A, const GeoShape* B);
 
-  //	Returns the volume of the shape, for mass inventory
+  //    Returns the volume of the shape, for mass inventory
   virtual double volume () const;
 
-  //	Returns the NOT shape type, as a string.
+  //    Returns the bonding box of the shape
+  virtual void extent (double& xmin, double& ymin, double& zmin,
+                       double& xmax, double& ymax, double& zmax) const;
+
+  //    Returns true if the shape contains the point, false otherwise
+  virtual bool contains (double x, double y, double z) const;
+
+  //    Returns the NOT shape type, as a string.
   virtual const std::string & type () const;
 
-  //	Returns the NOT shape type, as a coded integer.
+  //    Returns the NOT shape type, as a coded integer.
   virtual ShapeType typeID () const;
 
-  //	Returns the first operand in the subtraction
+  //    Returns the first operand in the subtraction
   const GeoShape* getOpA () const;
 
-  //	Returns the second operand in the subtraction
+  //    Returns the second operand in the subtraction
   const GeoShape* getOpB () const;
 
-  //	Executes a GeoShapeAction
+  //    Executes a GeoShapeAction
   virtual void exec (GeoShapeAction *action) const;
 
-  //	For type identification.
+  //    For type identification.
   static const std::string& getClassType ();
 
-  //	For type identification.
+  //    For type identification.
   static ShapeType getClassTypeID ();
 
  protected:
@@ -49,12 +57,15 @@ class GeoShapeSubtraction : public GeoShape
   GeoShapeSubtraction(const GeoShapeSubtraction &right);
   GeoShapeSubtraction & operator=(const GeoShapeSubtraction &right);
 
-  //	The shape operand in the Subtraction operation
+  //    The first shape operand in the Subtraction operation
   const GeoShape* m_opA;
 
-  //	The shape operand in the Subtraction operation
+  //    The second shape operand in the Subtraction operation
   const GeoShape* m_opB;
 
+  //    Cached volume
+  mutable double fVolume = -1.;
+
   static const std::string s_classType;
   static const ShapeType s_classTypeID;
 
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoShapeUnion.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoShapeUnion.h
index 0661ffcda6d0bf4011b73202cf0da2193d18132a..0fc3581fdf997d98f99ae4b199f1bf1e534209bf 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoShapeUnion.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoShapeUnion.h
@@ -14,30 +14,38 @@
 class GeoShapeUnion : public GeoShape
 {
  public:
+  //    Constructor taking two shape operands
   GeoShapeUnion (const GeoShape* A, const GeoShape* B);
 
-  //	Returns the volume of the shape, for mass inventory
+  //    Returns the volume of the shape, for mass inventory
   virtual double volume () const;
 
-  //	Returns the OR shape type, as a string.
+  //    Returns the bonding box of the shape
+  virtual void extent (double& xmin, double& ymin, double& zmin,
+                       double& xmax, double& ymax, double& zmax) const;
+
+  //    Returns true if the shape contains the point, false otherwise
+  virtual bool contains (double x, double y, double z) const;
+
+  //    Returns the OR shape type, as a string.
   virtual const std::string & type () const;
 
-  //	Returns the OR shape type, as a coded integer.
+  //    Returns the OR shape type, as a coded integer.
   virtual ShapeType typeID () const;
 
-  //	Returns the first operand being ORed
+  //    Returns the first operand being ORed
   const GeoShape* getOpA () const;
 
-  //	Returns the second operand being ORed.
+  //    Returns the second operand being ORed.
   const GeoShape* getOpB () const;
 
-  //	Executes a GeoShapeAction
+  //    Executes a GeoShapeAction
   virtual void exec (GeoShapeAction *action) const;
 
-  //	For type identification.
+  //    For type identification.
   static const std::string& getClassType ();
 
-  //	For type identification.
+  //    For type identification.
   static ShapeType getClassTypeID ();
 
  protected:
@@ -47,12 +55,15 @@ class GeoShapeUnion : public GeoShape
   GeoShapeUnion(const GeoShapeUnion &right);
   GeoShapeUnion & operator=(const GeoShapeUnion &right);
 
-  //	The first shape operand in the OR operation.
+  //    The first shape operand in the OR operation
   const GeoShape* m_opA;
 
-  //	The second shape operand in the OR operation.
+  //    The second shape operand in the OR operation
   const GeoShape* m_opB;
 
+  //    Cached volume
+  mutable double fVolume = -1.;
+
   static const std::string s_classType;
   static const ShapeType s_classTypeID;
 
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoSimplePolygonBrep.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoSimplePolygonBrep.h
index 288a47ddd6c4bcb072542484eced7c403ffcba65..74b00560a1749e43d2e8953974e98f3b5511669d 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoSimplePolygonBrep.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoSimplePolygonBrep.h
@@ -11,7 +11,7 @@
  * @brief This shape represents a BREP solid consisting
  *      of two equivalent polygonial faces perpendicular to Z axis.
  *      The polygones are described by array of (x,y) vertices,
- *      The solid is considered valid if the number of polygon vertices >=3 
+ *      The solid is considered valid if the number of polygon vertices >=3
  *
  *      Constructor parameter is a half length along Z axis
  */
@@ -19,30 +19,56 @@
 #include <vector>
 #include "GeoModelKernel/GeoShape.h"
 
-
 class GeoSimplePolygonBrep : public GeoShape
 {
  public:
+  //    Constructor for the BREP.  Note that the constructor
+  //    does not fully build this object. The BREP is not valid
+  //    until at least three vertices have been added.
   GeoSimplePolygonBrep(double dz);
 
+  //    Returns the volume of the shape, for mass inventory
   virtual double volume() const;
 
+  //    Returns the bonding box of the shape
+  virtual void extent (double& xmin, double& ymin, double& zmin,
+                       double& xmax, double& ymax, double& zmax) const;
+
+  //    Returns true if the shape contains the point, false otherwise.
+  virtual bool contains (double x, double y, double z) const;
+
+  //    Returns the BREP shape type, as a string
   virtual const std::string& type() const;
+
+  //    Returns the BREP shape type, as a coded integer
   virtual ShapeType typeID() const;
 
+  //    For type identification
+  static const std::string& getClassType();
+
+  //    For type identification
+  static ShapeType getClassTypeID();
+
+  //    Executes a GeoShapeAction
+  virtual void exec(GeoShapeAction *action) const;
+
+  //    Add another vertex to the polygon. A minimum of three
+  //    vertices are required to create a valid solid.
   void addVertex(double XVertex, double YVertex);
+
+  //    Returns the number of vertices in the polygon
   unsigned int getNVertices() const;
 
+  //    True if the polygon has at least three vertices, false otherwise
   bool isValid () const;
 
+  //    Returns X coordinate of the specified vertex
   const double & getXVertex(unsigned int i) const;
-  const double & getYVertex(unsigned int i) const;
-
-  virtual void exec(GeoShapeAction *action) const;
 
-  static const std::string& getClassType();
-  static ShapeType getClassTypeID();
+  //    Returns Y coordinate of the specified vertex
+  const double & getYVertex(unsigned int i) const;
 
+  //    Half-length along Z axis
   const double& getDZ() const;
 
  protected:
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTessellatedSolid.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTessellatedSolid.h
index 48640d6950ddda6ea7548c4bd1b4de52f0aded6a..c1d185ec5e7038aa7bcd1f9866135c43c2d58b56 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTessellatedSolid.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTessellatedSolid.h
@@ -11,24 +11,48 @@
 class GeoTessellatedSolid : public GeoShape
 {
  public:
+  //    Constructor for the TESSELLATED SOLID.  Note that the constructor
+  //    does not fully build this object.  The TESSELLATED SOLID is not valid
+  //    until at least four facets have been added.
   GeoTessellatedSolid();
 
+  //    Returns the volume of the shape, for mass inventory
   virtual double volume() const;
 
+  //    Returns the bonding box of the shape
+  virtual void extent (double& xmin, double& ymin, double& zmin,
+                       double& xmax, double& ymax, double& zmax) const;
+
+  //    Returns true if the shape contains the point, false otherwise.
+  virtual bool contains (double x, double y, double z) const;
+
+  //    Returns the TESSELLATED SOLID shape type, as a string
   virtual const std::string& type() const;
+
+  //    Returns the TESSELLATED SOLID shape type, as a coded integer
   virtual ShapeType typeID() const;
 
+  //    For type identification
   static const std::string& getClassType();
+
+  //    For type identification
   static ShapeType getClassTypeID();
 
+  //    Executes a GeoShapeAction
   virtual void exec(GeoShapeAction *action) const;
 
+  //    Add another facet to the tessellated solid. A minimum of four
+  //    facets are required to create a valid tesselated solid.
   void addFacet(GeoFacet*);
+
+  //    Returns specified facet
   GeoFacet* getFacet(size_t) const;
+
+  //    Returns the number of facets
   size_t getNumberOfFacets() const;
 
-  // True if the tessellated solid has at least four facets (tetrahedron).
-  // False otherwise.
+  // True if the tessellated solid has at least four facets (tetrahedron),
+  // false otherwise.
   bool isValid () const;
 
  protected:
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTorus.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTorus.h
index 54c4747d7258edf7af781ce440a8abf1a27d98d3..dce51eddbd5ddc714d2c790abae6f6e8f579a7c7 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTorus.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTorus.h
@@ -13,20 +13,44 @@ class GeoTorus : public GeoShape
  public:
   GeoTorus (double Rmin, double Rmax, double Rtor, double SPhi, double DPhi);
 
+  //    Returns the volume of the shape, for mass inventory.
   virtual double volume () const;
 
+  //    Returns the bonding box of the shape.
+  virtual void extent (double& xmin, double& ymin, double& zmin,
+                       double& xmax, double& ymax, double& zmax) const;
+
+  //    Returns true if the shape contains the point, false otherwise.
+  virtual bool contains (double x, double y, double z) const;
+
+  //    Returns the TORUS shape type, as a string.
   virtual const std::string & type () const;
+
+  //    Returns the TORUS 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 ();
-  static ShapeType getClassTypeID ();
 
-  virtual void exec (GeoShapeAction *action) const;
+  //    For type identification.
+  static ShapeType getClassTypeID ();
 
+  //    Inner radius.
   const double& getRMin () const;
+
+  //    Outer radius.
   const double& getRMax () const;
+
+  //    Radius of the torus.
   const double& getRTor () const;
+
+  //    Starting angle of the segment in radians.
   const double& getSPhi () const;
+
+  //    Delta angle of the segment in radians.
   const double& getDPhi () const;
 
  protected:
@@ -81,5 +105,4 @@ inline const double& GeoTorus::getDPhi () const
   return m_dPhi;
 }
 
-
 #endif
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTrap.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTrap.h
index 2f47b45bbed0140cfa7829d7be9dd12f640c5b20..f36693638d1fe43025569a7254ad04b76f5f3b44 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTrap.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTrap.h
@@ -18,55 +18,62 @@ class GeoTrap : public GeoShape
  public:
   GeoTrap (double ZHalfLength, double Theta, double Phi, double Dydzn, double Dxdyndzn, double Dxdypdzn, double Angleydzn, double Dydzp, double Dxdyndzp, double Dxdypdzp, double Angleydzp);
 
-  //	Returns the volume of the shape, for mass inventory
+  //    Returns the volume of the shape, for mass inventory.
   virtual double volume () const;
-  
-  //	Returns the TRAP shape type, as a string.
+
+  //    Returns the bonding box of the shape.
+  virtual void extent (double& xmin, double& ymin, double& zmin,
+                       double& xmax, double& ymax, double& zmax) const;
+
+  //    Returns true if the shape contains the point, false otherwise.
+  virtual bool contains (double x, double y, double z) const;
+
+  //    Returns the TRAP shape type, as a string.
   virtual const std::string & type () const;
-  
-  //	Returns the TRAP shape type, as a coded integer.
+
+  //    Returns the TRAP shape type, as a coded integer.
   virtual ShapeType typeID () const;
-  
-  //	Executes a GeoShapeAction
+
+  //    Executes a GeoShapeAction.
   virtual void exec (GeoShapeAction *action) const;
-  
-  //	For type identification.
+
+  //    For type identification.
   static const std::string& getClassType ();
-  
-  //	For type identification.
+
+  //    For type identification.
   static ShapeType getClassTypeID ();
-  
-  //	Z half length.
+
+  //    Z half length.
   const double& getZHalfLength () const;
-  
-  //	Polar (theta) angle.
+
+  //    Polar (theta) angle.
   const double& getTheta () const;
-  
-  //	Azimuthal (phi) angle.
+
+  //    Azimuthal (phi) angle.
   const double& getPhi () const;
-  
-  //	Y half length at -z.
+
+  //    Y half length at -z.
   const double& getDydzn () const;
-  
-  //	X half length at -z, -y.
+
+  //    X half length at -z, -y.
   const double& getDxdyndzn () const;
-  
-  //	X half length at -z, +y
+
+  //    X half length at -z, +y.
   const double& getDxdypdzn () const;
-  
+
   const double& getAngleydzn () const;
-  
-  //	Y half length at +z.
+
+  //    Y half length at +z.
   const double& getDydzp () const;
-  
-  //	X half length at +z, -y
+
+  //    X half length at +z, -y.
   const double& getDxdyndzp () const;
-  
-  //	X half length at +z, +y
+
+  //    X half length at +z, +y.
   const double& getDxdypdzp () const;
-  
+
   const double& getAngleydzp () const;
-  
+
  protected:
   virtual ~GeoTrap();
 
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTrd.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTrd.h
index 6bfa31c4ec152729652fdba21ea32bea9de1ebb4..bb28d149a6d80364f686d94868f84ad26ce03280 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTrd.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTrd.h
@@ -12,37 +12,44 @@ class GeoTrd : public GeoShape
  public:
   GeoTrd (double XHalfLength1, double XHalfLength2, double YHalfLength1, double YHalfLength2, double ZHalfLength);
 
-  //	Returns the volume of the shape, for mass inventory
+  //    Returns the volume of the shape, for mass inventory.
   virtual double volume () const;
-  
-  //	Returns the TRD shape type, as a string.
+
+  //    Returns the bonding box of the shape.
+  virtual void extent (double& xmin, double& ymin, double& zmin,
+                       double& xmax, double& ymax, double& zmax) const;
+
+  //    Returns true if the shape contains the point, false otherwise.
+  virtual bool contains (double x, double y, double z) const;
+
+  //    Returns the TRD shape type, as a string.
   virtual const std::string & type () const;
-  
-  //	Returns the TRD shape type, as a coded integer.
+
+  //    Returns the TRD shape type, as a coded integer.
   virtual ShapeType typeID () const;
-  
-  //	Executes a GeoShapeAction
+
+  //    Executes a GeoShapeAction.
   virtual void exec (GeoShapeAction *action) const;
-  
-  //	For type identification.
+
+  //    For type identification.
   static const std::string& getClassType ();
-  
-  //	For type identification.
+
+  //    For type identification.
   static ShapeType getClassTypeID ();
-  
-  //	Half length in the x-direction at -dz.
+
+  //    Half length in the x-direction at -dz.
   const double& getXHalfLength1 () const;
-  
-  //	Half length in the x-direction at +dz
+
+  //    Half length in the x-direction at +dz.
   const double& getXHalfLength2 () const;
-  
-  //	Half-length in the y direction at +dz.
+
+  //    Half-length in the y direction at +dz.
   const double& getYHalfLength1 () const;
-  
-  //	Half-length in the y direction at -dz
+
+  //    Half-length in the y direction at -dz.
   const double& getYHalfLength2 () const;
-  
-  //	Half-length in the z direction.
+
+  //    Half-length in the z direction.
   const double& getZHalfLength () const;
 
  protected:
@@ -51,7 +58,7 @@ class GeoTrd : public GeoShape
  private:
   GeoTrd(const GeoTrd &right);
   GeoTrd & operator=(const GeoTrd &right);
-  
+
   static const std::string s_classType;
   static const ShapeType s_classTypeID;
 
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTube.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTube.h
index 4354cfd4baa437c229ee25ac76cad45709231568..9cb7f9084edbe700128c6b8f2977c0a80d11c3c6 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTube.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTube.h
@@ -12,37 +12,44 @@ class GeoTube : public GeoShape
  public:
   GeoTube (double RMin, double RMax, double ZHalfLength);
 
-  //	Returns the volume of the shape, for mass inventory
+  //    Returns the volume of the shape, for mass inventory.
   virtual double volume () const;
-  
-  //	Returns the TUBE shape type, as a string.
+
+  //    Returns the bonding box of the shape.
+  virtual void extent (double& xmin, double& ymin, double& zmin,
+                       double& xmax, double& ymax, double& zmax) const;
+
+  //    Returns true if the shape contains the point, false otherwise.
+  virtual bool contains (double x, double y, double z) const;
+
+  //    Returns the TUBE shape type, as a string.
   virtual const std::string & type () const;
-  
-  //	Returns the TUBE shape type, as a coded integer.
+
+  //    Returns the TUBE shape type, as a coded integer.
   virtual ShapeType typeID () const;
-  
-  //	Executes a GeoShapeAction
+
+  //    Executes a GeoShapeAction.
   virtual void exec (GeoShapeAction *action) const;
-  
-  //	For type identification.
+
+  //    For type identification.
   static const std::string& getClassType ();
-  
-  //	For type identification.
+
+  //    For type identification.
   static ShapeType getClassTypeID ();
-  
-  //	Minimum (inner) tube radius.
+
+  //    Minimum (inner) tube radius.
   const double& getRMin () const;
-  
-  //	Maximum (outer) tube radius.
+
+  //    Maximum (outer) tube radius.
   const double& getRMax () const;
-  
-  //	Tube half-length in the z direction.
+
+  //    Tube half-length in the z direction.
   const double& getZHalfLength () const;
-  
+
  protected:
   //## Destructor (generated)
   virtual ~GeoTube();
-  
+
  private:
   GeoTube(const GeoTube &right);
   GeoTube & operator=(const GeoTube &right);
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTubs.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTubs.h
index 68f8399fb5b2ed27cee402f3d3fc8520238cf1d6..7aab48ed526ad3be0d20544b0415d9b3d5e49607 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTubs.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTubs.h
@@ -12,37 +12,44 @@ class GeoTubs : public GeoShape
  public:
   GeoTubs (double RMin, double RMax, double ZHalfLength, double SPhi, double DPhi);
 
-  //	Returns the volume of the shape, for mass inventory
+  //    Returns the volume of the shape, for mass inventory.
   virtual double volume () const;
 
-  //	Returns the TUBS shape type, as a string.
+  //    Returns the bonding box of the shape.
+  virtual void extent (double& xmin, double& ymin, double& zmin,
+                       double& xmax, double& ymax, double& zmax) const;
+
+  //    Returns true if the shape contains the point, false otherwise.
+  virtual bool contains (double x, double y, double z) const;
+
+  //    Returns the TUBS shape type, as a string.
   virtual const std::string & type () const;
 
-  //	Returns the TUBS shape type, as a coded integer.
+  //    Returns the TUBS shape type, as a coded integer.
   virtual ShapeType typeID () const;
 
-  //	Executes a GeoShapeAction
+  //    Executes a GeoShapeAction.
   virtual void exec (GeoShapeAction *action) const;
 
-  //	For type identification.
+  //    For type identification.
   static const std::string& getClassType ();
 
-  //	For type identification.
+  //    For type identification.
   static ShapeType getClassTypeID ();
 
-  //	Minimum (inner) tube section radius.
+  //    Minimum (inner) tube section radius.
   const double& getRMin () const;
 
-  //	Maximum (outer) tube section radius.
+  //    Maximum (outer) tube section radius.
   const double& getRMax () const;
 
-  //	Tube section half-length in the z direction.
+  //    Tube section half-length in the z direction.
   const double& getZHalfLength () const;
 
-  //	Starting angle of the tube section in radians.
+  //    Starting angle of the tube section in radians.
   const double& getSPhi () const;
 
-  //	Delta angle of the tube section in radians.
+  //    Delta angle of the tube section in radians.
   const double& getDPhi () const;
 
  protected:
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTwistedTrap.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTwistedTrap.h
index e88cc4665ab9bea7e634b84ee4e96ccb47ede466..a7e6b5261bf767a4a5acc85b6d2a37a00ed820a9 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTwistedTrap.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTwistedTrap.h
@@ -43,13 +43,13 @@
 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
@@ -63,25 +63,32 @@ class GeoTwistedTrap : public GeoShape
                  double  pAlph    // tilt angle
   );
 
-  //	Returns the volume of the shape, for mass inventory
+  //    Returns the volume of the shape, for mass inventory
   virtual double volume () const;
-  
-  //	Returns the TWISTED TRAP shape type, as a string.
+
+  //    Returns the bonding box of the shape.
+  virtual void extent (double& xmin, double& ymin, double& zmin,
+                       double& xmax, double& ymax, double& zmax) const;
+
+  //    Returns true if the shape contains the point, false otherwise.
+  virtual bool contains (double x, double y, double z) 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.
+
+  //    Returns the TWISTED TRAP shape type, as a coded integer.
   virtual ShapeType typeID () const;
-  
-  //	Executes a GeoShapeAction
+
+  //    Executes a GeoShapeAction
   virtual void exec (GeoShapeAction *action) const;
-  
-  //	For type identification.
+
+  //    For type identification.
   static const std::string& getClassType ();
-  
-  //	For type identification.
+
+  //    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 ; }
@@ -93,7 +100,7 @@ class GeoTwistedTrap : public GeoShape
   inline double getTheta()        const { return m_theta ; }
   inline double getPhi()          const { return m_phi ; }
   inline double getTiltAngleAlpha() const { return m_alph ; }
-  
+
  protected:
   virtual ~GeoTwistedTrap();
 
@@ -115,7 +122,7 @@ class GeoTwistedTrap : public GeoShape
   double m_dz;        // Half-length along the z axis
   double m_alph ;
   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;
@@ -131,4 +138,5 @@ inline ShapeType GeoTwistedTrap::getClassTypeID ()
 {
   return s_classTypeID;
 }
+
 #endif
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoUnidentifiedShape.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoUnidentifiedShape.h
index 22428ec6a60b3ba189d8d73f905c1c4da247b6ce..8f021aee87077af65fcf6b5bd3f4b7e9b3bd2282 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoUnidentifiedShape.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoUnidentifiedShape.h
@@ -8,7 +8,7 @@
 // This class is essentially a bundle of data with no fixed format.
 // ASCII information can be added and retreived. It serves as a
 // proxy for shapes that are not part of the kernel, i.e. extender-
-// class shapes. 
+// class shapes.
 //
 // Joe Boudreau 2020
 //
@@ -19,17 +19,17 @@
 #include "GeoModelKernel/GeoShape.h"
 #include "GeoModelKernel/Query.h"
 #include <string>
+
 class GeoUnidentifiedShape: public GeoShape {
 
  public:
-
   // Constructor:
   GeoUnidentifiedShape(const std::string & name);
-  
+
   // Constructor:
   GeoUnidentifiedShape(const std::string & name, const std::string & asciiData);
-  
-  // Constructor with volume: 
+
+  // Constructor with volume:
   GeoUnidentifiedShape(const std::string & name, const std::string & asciiData, double volume);
 
   // Returns the user-provided name of the volume (eg "MySpecialShape");
@@ -37,32 +37,39 @@ class GeoUnidentifiedShape: public GeoShape {
 
   // Returns the ascii data associated with this object (possibly empty);
   const std::string & asciiData() const;
-  
+
   // Returns the volume of the shape, for mass inventory
   virtual double volume () const;
-  
+
+  // Returns the bonding box of the shape
+  virtual void extent (double& xmin, double& ymin, double& zmin,
+                       double& xmax, double& ymax, double& zmax) const;
+
+  // Returns true if the shape contains the point, false otherwise
+  virtual bool contains (double x, double y, double z) const;
+
   // Returns the shape type, as a string.
   virtual const std::string & type () const;
-  
+
   // Returns the shape type, as an coded integer.
   virtual ShapeType typeID () const;
 
-  //	For type identification.
+  // For type identification.
   static const std::string& getClassType ();
 
-  //	For type identification.,
+  // For type identification.,
   static ShapeType getClassTypeID ();
-  
-  
+
+
   // Executes a GeoShapeAction
   virtual void exec (GeoShapeAction *action) const;
 
 
-  
+
  protected:
   // Destructor:
   ~GeoUnidentifiedShape();
-  
+
  private:
 
   const std::string         _name;
@@ -70,7 +77,7 @@ class GeoUnidentifiedShape: public GeoShape {
   const Query<double>       _volume;
   static const std::string         _classType;
   static const ShapeType           _classTypeID;
-  
+
 };
 
 inline const std::string& GeoUnidentifiedShape::getClassType ()
@@ -82,4 +89,5 @@ inline ShapeType GeoUnidentifiedShape::getClassTypeID ()
 {
   return _classTypeID;
 }
+
 #endif
diff --git a/GeoModelCore/GeoModelKernel/src/GeoBox.cxx b/GeoModelCore/GeoModelKernel/src/GeoBox.cxx
index de07cf79cc71a3431b265ece87e72d42d35c4fcd..7c881ed60e628d2c730d0fc8a412bacaeed72464 100755
--- a/GeoModelCore/GeoModelKernel/src/GeoBox.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoBox.cxx
@@ -24,6 +24,25 @@ double GeoBox::volume () const
   return 8.0 * m_xHalfLength * m_yHalfLength * m_zHalfLength;
 }
 
+void GeoBox::extent (double& xmin, double& ymin, double& zmin,
+                     double& xmax, double& ymax, double& zmax) const
+{
+  xmin = -m_xHalfLength;
+  ymin = -m_yHalfLength;
+  zmin = -m_zHalfLength;
+  xmax = m_xHalfLength;
+  ymax = m_yHalfLength;
+  zmax = m_zHalfLength;
+}
+
+bool GeoBox::contains (double x, double y, double z) const
+{
+  double distx = std::abs(x) - m_xHalfLength;
+  double disty = std::abs(y) - m_yHalfLength;
+  double distz = std::abs(z) - m_zHalfLength;
+  return (std::max(std::max(distx, disty), distz) <= 0.0);
+}
+
 const std::string & GeoBox::type () const
 {
   return s_classType;
diff --git a/GeoModelCore/GeoModelKernel/src/GeoCons.cxx b/GeoModelCore/GeoModelKernel/src/GeoCons.cxx
index 1082bfb271db4c122e57bf9ccd78e320976d51d4..10d3959a9430d767d8c454d41b0428ec407b7e28 100755
--- a/GeoModelCore/GeoModelKernel/src/GeoCons.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoCons.cxx
@@ -4,6 +4,7 @@
 
 #include "GeoModelKernel/GeoCons.h"
 #include "GeoModelKernel/GeoShapeAction.h"
+#include <cmath>
 
 const std::string GeoCons::s_classType = "Cons";
 const ShapeType GeoCons::s_classTypeID = 0x11;
@@ -25,9 +26,40 @@ GeoCons::~GeoCons()
 
 double GeoCons::volume () const
 {
-  return (m_dPhi * (1./3)) * m_dZ * (m_rMax1 * m_rMax1 + m_rMax2 * m_rMax2 + m_rMax1 * m_rMax2
-                                 - m_rMin1 * m_rMin1 - m_rMin2 * m_rMin2 -
-                                 m_rMin1 * m_rMin2);
+  return (m_dPhi * (1./3.)) * m_dZ *
+    (m_rMax1 * m_rMax1 + m_rMax2 * m_rMax2 + m_rMax1 * m_rMax2 -
+     m_rMin1 * m_rMin1 - m_rMin2 * m_rMin2 - m_rMin1 * m_rMin2);
+}
+
+void GeoCons::extent (double& xmin, double& ymin, double& zmin,
+                      double& xmax, double& ymax, double& zmax) const
+{
+  double rmin = std::min(m_rMin1, m_rMin2);
+  double rmax = std::max(m_rMax1, m_rMax2);
+  GeoShape::diskExtent(rmin, rmax, m_sPhi, m_dPhi, xmin, ymin, xmax, ymax);
+  zmin =-m_dZ;
+  zmax = m_dZ;
+}
+
+bool GeoCons::contains (double x, double y, double z) const
+{
+#ifndef M_PI
+  constexpr double M_PI = 3.14159265358979323846;
+#endif
+  if (std::abs(z) - m_dZ > 0.0) return false;
+  double t = 0.5 * (1.0 + z / m_dZ);
+  double rmin = m_rMin1 + (m_rMin2 - m_rMin1) * t;
+  double rmax = m_rMax1 + (m_rMax2 - m_rMax1) * t;
+  double rr = x * x + y * y;
+  if (rr > rmax * rmax || rr < rmin * rmin) return false;
+  if (m_dPhi >= 2.0 * M_PI) return true;
+
+  GeoTrf::Vector2D r(x, y);
+  GeoTrf::Vector2D ns(std::sin(m_sPhi), -std::cos(m_sPhi));
+  GeoTrf::Vector2D ne(-std::sin(m_sPhi + m_dPhi), std::cos(m_sPhi + m_dPhi));
+  double ds = ns.dot(r);
+  double de = ne.dot(r);
+  return (m_dPhi <= M_PI) ? (ds <= 0 && de <= 0) : (ds <= 0 || de <= 0);
 }
 
 const std::string & GeoCons::type () const
diff --git a/GeoModelCore/GeoModelKernel/src/GeoEllipticalTube.cxx b/GeoModelCore/GeoModelKernel/src/GeoEllipticalTube.cxx
index d0848e1111b67561cdefd7ad37205c3c2429c5fe..91b52fc9e2e0a3d82952a3c6a3c9465b92000435 100755
--- a/GeoModelCore/GeoModelKernel/src/GeoEllipticalTube.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoEllipticalTube.cxx
@@ -2,15 +2,14 @@
   Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
 */
 
-#include <cmath>
 #include "GeoModelKernel/GeoPolyhedrizeAction.h"
 #include "GeoModelKernel/GeoShapeAction.h"
+#include <cmath>
 
 #include "GeoModelKernel/GeoEllipticalTube.h"
 const std::string GeoEllipticalTube::s_classType = "EllipticalTube";
 const ShapeType GeoEllipticalTube::s_classTypeID = 0x22;
 
-
 GeoEllipticalTube::GeoEllipticalTube(double XHalfLength, double YHalfLength, double ZHalfLength)
   : m_xHalfLength(XHalfLength)
   , m_yHalfLength(YHalfLength)
@@ -18,7 +17,6 @@ GeoEllipticalTube::GeoEllipticalTube(double XHalfLength, double YHalfLength, dou
 {
 }
 
-
 GeoEllipticalTube::~GeoEllipticalTube()
 {
 }
@@ -26,9 +24,27 @@ GeoEllipticalTube::~GeoEllipticalTube()
 double GeoEllipticalTube::volume () const
 {
 #ifndef M_PI
-  double M_PI = acos (-1.0);
+  constexpr double M_PI = 3.14159265358979323846;
 #endif
-  return 2 * M_PI * m_xHalfLength * m_yHalfLength * m_zHalfLength;
+  return 2.0 * M_PI * m_xHalfLength * m_yHalfLength * m_zHalfLength;
+}
+
+void GeoEllipticalTube::extent (double& xmin, double& ymin, double& zmin,
+                                double& xmax, double& ymax, double& zmax) const
+{
+  xmin =-m_xHalfLength;
+  ymin =-m_yHalfLength;
+  zmin =-m_zHalfLength;
+  xmax = m_xHalfLength;
+  ymax = m_yHalfLength;
+  zmax = m_zHalfLength;
+}
+
+bool GeoEllipticalTube::contains (double x, double y, double z) const
+{
+  if (std::abs(z) - m_zHalfLength > 0.0) return false;
+  return ((x * x) / (m_xHalfLength * m_xHalfLength) +
+          (y * y) / (m_yHalfLength * m_yHalfLength) <= 1.0);
 }
 
 const std::string & GeoEllipticalTube::type () const
@@ -45,4 +61,3 @@ void GeoEllipticalTube::exec (GeoShapeAction *action) const
 {
   action->handleEllipticalTube(this);
 }
-
diff --git a/GeoModelCore/GeoModelKernel/src/GeoGenericTrap.cxx b/GeoModelCore/GeoModelKernel/src/GeoGenericTrap.cxx
index 346d99f6a15f6b081bd604e46d507839d7a666d5..a347dae785acf14fb1d2df3c70dd34f2a8ad2f03 100755
--- a/GeoModelCore/GeoModelKernel/src/GeoGenericTrap.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoGenericTrap.cxx
@@ -21,18 +21,57 @@ GeoGenericTrap::~GeoGenericTrap()
 double GeoGenericTrap::volume() const
 {
   // diagonals
-  GeoTwoVector A = m_vertices[3] - m_vertices[1];
-  GeoTwoVector B = m_vertices[2] - m_vertices[0];
-  GeoTwoVector C = m_vertices[7] - m_vertices[5];
-  GeoTwoVector D = m_vertices[6] - m_vertices[4];
+  GeoTwoVector a = m_vertices[3] - m_vertices[1];
+  GeoTwoVector b = m_vertices[2] - m_vertices[0];
+  GeoTwoVector c = m_vertices[7] - m_vertices[5];
+  GeoTwoVector d = m_vertices[6] - m_vertices[4];
 
   // kross products
-  double AB = A.x()*B.y() - A.y()*B.x();
-  double CD = C.x()*D.y() - C.y()*D.x();
-  double AD = A.x()*D.y() - A.y()*D.x();
-  double CB = C.x()*B.y() - C.y()*B.x();
+  double ab = a.x()*b.y() - a.y()*b.x();
+  double cd = c.x()*d.y() - c.y()*d.x();
+  double ad = a.x()*d.y() - a.y()*d.x();
+  double cb = c.x()*b.y() - c.y()*b.x();
 
-  return m_zHalfLength * ((AB + CD)* (1./3.) + (AD + CB)* (1./6.));
+  return m_zHalfLength * ((ab + cd)* (1./3.) + (ad + cb)* (1./6.));
+}
+
+void GeoGenericTrap::extent (double& xmin, double& ymin, double& zmin,
+                             double& xmax, double& ymax, double& zmax) const
+{
+  xmin = xmax = m_vertices[0].x();
+  ymin = ymax = m_vertices[0].y();
+  int nv = m_vertices.size();
+  for (int i = 1; i < nv; ++i)
+  {
+    if (xmin > m_vertices[i].x()) xmin = m_vertices[i].x();
+    if (xmax < m_vertices[i].x()) xmax = m_vertices[i].x();
+    if (ymin > m_vertices[i].y()) ymin = m_vertices[i].y();
+    if (ymax < m_vertices[i].y()) ymax = m_vertices[i].y();
+  }
+  zmin =-m_zHalfLength;
+  zmax = m_zHalfLength;
+}
+
+bool GeoGenericTrap::contains (double x, double y, double z) const
+{
+  if (std::abs(z) - m_zHalfLength > 0.0) return false;
+  double t = 0.5 * (1.0 + z / m_zHalfLength);
+  GeoTrf::Vector2D v[4];
+  v[0] = m_vertices[0] + (m_vertices[4] - m_vertices[0]) * t;
+  v[1] = m_vertices[1] + (m_vertices[5] - m_vertices[1]) * t;
+  v[2] = m_vertices[2] + (m_vertices[6] - m_vertices[2]) * t;
+  v[3] = m_vertices[3] + (m_vertices[7] - m_vertices[3]) * t;
+  int nv = 4;
+  bool in = false;
+  for (int i = 0, k = nv - 1; i < nv; k = i++)
+  {
+    if ((v[i].y() > y) != (v[k].y() > y))
+    {
+      double ctg = (v[k].x() - v[i].x()) / (v[k].y() - v[i].y());
+      in ^= (x < (y - v[i].y()) * ctg + v[i].x());
+    }
+  }
+  return in;
 }
 
 const std::string& GeoGenericTrap::type() const
@@ -59,4 +98,3 @@ const GeoGenericTrapVertices& GeoGenericTrap::getVertices() const
 {
   return m_vertices;
 }
-
diff --git a/GeoModelCore/GeoModelKernel/src/GeoPara.cxx b/GeoModelCore/GeoModelKernel/src/GeoPara.cxx
index fbda7108c1c4e8b8b7f1ebcc51df66f00073c709..d515352630c9cb37fea3f9256052d730c90d638e 100755
--- a/GeoModelCore/GeoModelKernel/src/GeoPara.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoPara.cxx
@@ -27,6 +27,42 @@ double GeoPara::volume () const
   return 8.0 * m_xHalfLength * m_yHalfLength * m_zHalfLength;
 }
 
+void GeoPara::extent (double& xmin, double& ymin, double& zmin,
+                      double& xmax, double& ymax, double& zmax) const
+{
+  double dx = m_xHalfLength;
+  double dy = m_yHalfLength;
+  double dz = m_zHalfLength;
+
+  double x0 = dz * std::tan(m_theta) * std::cos(m_phi);
+  double x1 = dy * std::tan(m_alpha);
+  xmin = std::min(std::min(std::min(-x0-x1-dx, -x0+x1-dx), x0-x1-dx), x0+x1-dx);
+  xmax = std::max(std::max(std::max(-x0-x1+dx, -x0+x1+dx), x0-x1+dx), x0+x1+dx);
+
+  double y0 = dz * std::tan(m_theta) * std::sin(m_phi);
+  ymin = std::min(-y0-dy, y0-dy);
+  ymax = std::max(-y0+dy, y0+dy);
+
+  zmin = -dz;
+  zmax = dz;
+}
+
+bool GeoPara::contains (double x, double y, double z) const
+{
+  double cosPhi = std::cos(m_phi);
+  double sinPhi = std::sin(m_phi);
+  double tanTheta = std::tan(m_theta);
+  double tanAlpha = std::tan(m_alpha);
+
+  double z0 = z;
+  double y0 = y - z0 * tanTheta * sinPhi;
+  double x0 = x - y0 * tanAlpha - z0 * tanTheta * cosPhi;
+  double distx = std::abs(x0) - m_xHalfLength;
+  double disty = std::abs(y0) - m_yHalfLength;
+  double distz = std::abs(z0) - m_zHalfLength;
+  return (std::max(std::max(distx, disty), distz) <= 0.0);
+}
+
 const std::string & GeoPara::type () const
 {
   return s_classType;
@@ -39,6 +75,5 @@ ShapeType GeoPara::typeID () const
 
 void GeoPara::exec (GeoShapeAction *action) const
 {
-	action->handlePara(this);
+  action->handlePara(this);
 }
-
diff --git a/GeoModelCore/GeoModelKernel/src/GeoPcon.cxx b/GeoModelCore/GeoModelKernel/src/GeoPcon.cxx
index a7407b86c71b0855b93d08a1e94c330a6a52bca0..e7ad137ea9c82e8329f636ac92d01148963f86cd 100755
--- a/GeoModelCore/GeoModelKernel/src/GeoPcon.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoPcon.cxx
@@ -24,20 +24,75 @@ double GeoPcon::volume () const
 {
   if (!isValid ())
     throw std::runtime_error ("Volume requested for incomplete polycone");
-  double v = 0;
-  for (size_t s = 0; s < getNPlanes () - 1; s++)
-    {
-      double fDz = fabs (getZPlane (s + 1) - getZPlane (s)) / 2.0;
-      double fRmin1 = getRMinPlane (s + 1);
-      double fRmin2 = getRMinPlane (s);
-      double fRmax1 = getRMaxPlane (s + 1);
-      double fRmax2 = getRMaxPlane (s);
-      v +=
-	(m_dPhi * (1./3)) * fDz * (fRmax1 * fRmax1 + fRmax2 * fRmax2 +
-                                   fRmax1 * fRmax2 - fRmin1 * fRmin1 -
-                                   fRmin2 * fRmin2 - fRmin1 * fRmin2);
-    }
-  return v;
+  double vol = 0;
+  for (size_t k = 0; k < getNPlanes() - 1; ++k)
+  {
+    double fDz = std::abs(getZPlane(k + 1) - getZPlane(k)) / 2.0;
+    double fRmin1 = getRMinPlane(k + 1);
+    double fRmin2 = getRMinPlane(k);
+    double fRmax1 = getRMaxPlane(k + 1);
+    double fRmax2 = getRMaxPlane(k);
+    vol += (m_dPhi * (1./3)) * fDz *
+      (fRmax1 * fRmax1 + fRmax2 * fRmax2 + fRmax1 * fRmax2 -
+       fRmin1 * fRmin1 - fRmin2 * fRmin2 - fRmin1 * fRmin2);
+  }
+  return vol;
+}
+
+void GeoPcon::extent (double& xmin, double& ymin, double& zmin,
+                      double& xmax, double& ymax, double& zmax) const
+{
+  if (!isValid ())
+    throw std::runtime_error ("Extent requested for incomplete polycone");
+  double rmin = getRMinPlane(0);
+  double rmax = getRMaxPlane(0);
+  zmin = zmax = getZPlane(0);
+  for (size_t i = 1; i < getNPlanes(); ++i)
+  {
+    rmin = std::min(rmin, getRMinPlane(i));
+    rmax = std::max(rmax, getRMaxPlane(i));
+    zmin = std::min(zmin, getZPlane(i));
+    zmax = std::max(zmax, getZPlane(i));
+  }
+  GeoShape::diskExtent(rmin, rmax, getSPhi(), getDPhi(), xmin, ymin, xmax, ymax);
+}
+
+bool GeoPcon::contains (double x, double y, double z) const
+{
+#ifndef M_PI
+  constexpr double M_PI = 3.14159265358979323846;
+#endif
+  if (!isValid ()) return false;
+  size_t nz = getNPlanes();
+  if (z < getZPlane(0) || z > getZPlane(nz - 1)) return false;
+
+  if (m_dPhi < 2.0 * M_PI)
+  {
+    GeoTrf::Vector2D r(x, y);
+    GeoTrf::Vector2D ns(std::sin(m_sPhi), -std::cos(m_sPhi));
+    GeoTrf::Vector2D ne(-std::sin(m_sPhi + m_dPhi), std::cos(m_sPhi + m_dPhi));
+    double ds = ns.dot(r);
+    double de = ne.dot(r);
+    bool inwedge = (m_dPhi <= M_PI) ? (ds <= 0 && de <= 0) : (ds <= 0 || de <= 0);
+    if (!inwedge) return false;
+  }
+
+  double rr = x * x + y * y;
+  for (size_t k = 0; k < nz - 1; ++k)
+  {
+    double zmin = getZPlane(k);
+    double zmax = getZPlane(k + 1);
+    if (z < zmin || z > zmax || zmin == zmax) continue;
+    double t = (z - zmin) / (zmax - zmin);
+    double rmin1 = getRMinPlane(k);
+    double rmin2 = getRMinPlane(k + 1);
+    double rmin = rmin1 + (rmin2 - rmin1) * t;
+    double rmax1 = getRMaxPlane(k);
+    double rmax2 = getRMaxPlane(k + 1);
+    double rmax = rmax1 + (rmax2 - rmax1) * t;
+    if (rr <= rmax * rmax && rr >= rmin * rmin) return true;
+  }
+  return false;
 }
 
 const std::string & GeoPcon::type () const
@@ -61,4 +116,3 @@ void GeoPcon::exec (GeoShapeAction *action) const
 {
   action->handlePcon(this);
 }
-
diff --git a/GeoModelCore/GeoModelKernel/src/GeoPgon.cxx b/GeoModelCore/GeoModelKernel/src/GeoPgon.cxx
index 481740d865985a077396c2ace1ad61aaca519c3e..4412042ae2d4680baa87bb108a3e65daf7b225cb 100755
--- a/GeoModelCore/GeoModelKernel/src/GeoPgon.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoPgon.cxx
@@ -23,9 +23,12 @@ GeoPgon::~GeoPgon()
 
 double GeoPgon::volume () const
 {
+#ifndef M_PI
+  constexpr double M_PI = 3.14159265358979323846;
+#endif
   if (!isValid ())
     throw std::runtime_error ("Volume requested for incomplete polygon");
-  double v = 0.0;
+  double vol = 0.0;
   for (size_t k = 0; k < getNPlanes() - 1; ++k) {
     double z1 = getZPlane(k);
     double z2 = getZPlane(k + 1);
@@ -33,13 +36,123 @@ double GeoPgon::volume () const
     double a2 = getRMinPlane(k + 1);
     double b1 = getRMaxPlane(k);
     double b2 = getRMaxPlane(k + 1);
-    v += (b1*b1 + b1*b2 + b2*b2 - a1*a1 - a1*a2 - a2*a2) * (z2 - z1);
+    vol += (b1*b1 + b1*b2 + b2*b2 - a1*a1 - a1*a2 - a2*a2) * (z2 - z1);
   }
   int nsides = getNSides();
-  double alpha = m_dPhi/nsides;
-  double sinAlpha = sin(alpha);
-  double cosHalfAlpha = cos(0.5 * alpha);
-  return fabs(v) * nsides * sinAlpha / (cosHalfAlpha * cosHalfAlpha * 6.0);
+  double dphi = getDPhi();
+  if (dphi > 2.0 * M_PI) dphi = 2.0 * M_PI;
+  double alpha = dphi / nsides;
+  double sinAlpha = std::sin(alpha);
+  double cosHalfAlpha = std::cos(0.5 * alpha);
+  return std::abs(vol) * nsides * sinAlpha / (cosHalfAlpha * cosHalfAlpha * 6.0);
+}
+
+void GeoPgon::extent (double& xmin, double& ymin, double& zmin,
+                      double& xmax, double& ymax, double& zmax) const
+{
+#ifndef M_PI
+  constexpr double M_PI = 3.14159265358979323846;
+#endif
+  if (!isValid ())
+    throw std::runtime_error ("Extent requested for incomplete polygon");
+  double rmin = getRMinPlane(0);
+  double rmax = getRMaxPlane(0);
+  zmin = zmax = getZPlane(0);
+  for (size_t k = 1; k < getNPlanes(); ++k)
+  {
+    rmin = std::min(rmin, getRMinPlane(k));
+    rmax = std::max(rmax, getRMaxPlane(k));
+    zmin = std::min(zmin, getZPlane(k));
+    zmax = std::max(zmax, getZPlane(k));
+  }
+  size_t nsides = getNSides();
+  double dphi = getDPhi();
+  if (dphi >  2.0 * M_PI) dphi = 2.0 * M_PI;
+  if (dphi == 2.0 * M_PI) rmin = 0.0;
+  double alpha = dphi / nsides;
+  double invCosHalfAlpha = 1. / std::cos(0.5 * alpha);
+  rmin *= invCosHalfAlpha;
+  rmax *= invCosHalfAlpha;
+  double phi = getSPhi();
+  xmin = xmax = rmin * std::cos(phi);
+  ymin = ymax = rmin * std::sin(phi);
+  for (size_t k = 0; k <= nsides; ++k)
+  {
+    double x = rmax * std::cos(phi);
+    xmin = std::min(xmin, x);
+    xmax = std::max(xmax, x);
+    double y = rmax * std::sin(phi);
+    ymin = std::min(ymin, y);
+    ymax = std::max(ymax, y);
+    if (rmin > 0.)
+    {
+      double xx = rmin * std::cos(phi);
+      xmin = std::min(xmin, xx);
+      xmax = std::max(xmax, xx);
+      double yy = rmin * std::sin(phi);
+      ymin = std::min(ymin, yy);
+      ymax = std::max(ymax, yy);
+    }
+    phi += alpha;
+  }
+}
+
+bool GeoPgon::contains (double x, double y, double z) const
+{
+  if (!isValid ()) return false;
+  size_t nz = getNPlanes();
+  if (z < getZPlane(0) || z > getZPlane(nz - 1)) return false;
+
+  size_t nsides = getNSides();
+  double sphi = getSPhi();
+  double dphi = getDPhi();
+
+  double dangle = dphi / nsides; // sector angle 
+  double dhalfangle = 0.5 * dangle;
+  double cosHalfAngle = std::cos(dhalfangle);
+  double sinHalfAngle = std::sin(dhalfangle);
+  double cosAngle = (cosHalfAngle + sinHalfAngle) * (cosHalfAngle - sinHalfAngle);
+  double sinAngle = 2.0 * cosHalfAngle * sinHalfAngle;
+
+  double r = 0.0;
+  double rot = -(sphi + dhalfangle); // initial rotation
+  double cosRot = std::cos(rot);
+  double sinRot = std::sin(rot);
+  bool inwedge = false;
+  for (size_t iside = 0; iside < nsides; ++iside)
+  {
+    double xrot = x * cosRot - y * sinRot;
+    double yrot = x * sinRot + y * cosRot;
+    double dista = sinHalfAngle * xrot + cosHalfAngle * yrot;
+    double distb = sinHalfAngle * xrot - cosHalfAngle * yrot;
+    if (dista >= 0.0 && distb >= 0.0)
+    {
+      r = xrot;
+      inwedge = true;
+      break;
+    }
+    double cosTmp = cosRot;
+    double sinTmp = sinRot;
+    cosRot = cosTmp * cosAngle + sinTmp * sinAngle;
+    sinRot = sinTmp * cosAngle - cosTmp * sinAngle;
+  }
+  if (!inwedge) return false;
+
+  for (size_t k = 0; k < nz - 1; ++k)
+  {
+    double zmin = getZPlane(k);
+    double zmax = getZPlane(k + 1);
+    if (z < zmin || z > zmax || zmin == zmax) continue;
+    double t = (z - zmin) / (zmax - zmin);
+    double rmin1 = getRMinPlane(k);
+    double rmin2 = getRMinPlane(k + 1);
+    double rmin = rmin1 + (rmin2 - rmin1) * t;
+    double rmax1 = getRMaxPlane(k);
+    double rmax2 = getRMaxPlane(k + 1);
+    double rmax = rmax1 + (rmax2 - rmax1) * t;
+    if (r <= rmax && r >= rmin) return true;
+  }
+  return false;
 }
 
 const std::string & GeoPgon::type () const
diff --git a/GeoModelCore/GeoModelKernel/src/GeoShape.cxx b/GeoModelCore/GeoModelKernel/src/GeoShape.cxx
index 8142e0a0fb5fb1f87b4058fbc44f6518402446fc..96bd8c2002871ab03297ea3724071ee59abb6686 100755
--- a/GeoModelCore/GeoModelKernel/src/GeoShape.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoShape.cxx
@@ -7,6 +7,8 @@
 #include "GeoModelKernel/GeoShapeUnion.h"
 #include "GeoModelKernel/GeoShapeSubtraction.h"
 #include "GeoModelKernel/GeoShapeShift.h"
+#include <cmath>
+#include <cstdint>
 
 GeoShape::GeoShape ()
 {
@@ -16,6 +18,54 @@ GeoShape::~GeoShape()
 {
 }
 
+double GeoShape::volume () const
+{
+  constexpr int npoints = 1000000;     // number of random points
+  constexpr double expansion = 0.001;  // bounding box expansion
+  constexpr double f = 1./4294967296.; // 2^-32 - int to double conversion
+
+  // set up bonding box
+  double xmin = 0, ymin = 0, zmin = 0, xmax = 0, ymax = 0, zmax = 0;
+  extent(xmin, ymin, zmin, xmax, ymax, zmax);
+  double delx = (xmax - xmin) * expansion;
+  double dely = (ymax - ymin) * expansion;
+  double delz = (zmax - zmin) * expansion;
+  xmin -= delx;
+  ymin -= dely;
+  zmin -= delz;
+  xmax += delx;
+  ymax += dely;
+  zmax += delz;
+
+  uint32_t y = 2463534242; // seed for random number generation
+  int icount = 0; // counter of inside points
+  for (auto i = 0; i < npoints; ++i)
+  {
+    // generate three random numbers
+    uint32_t x = y;
+    x ^= x << 13;
+    x ^= x >> 17;
+    x ^= x << 5;
+    double randx = x * f;
+    x ^= x << 13;
+    x ^= x >> 17;
+    x ^= x << 5;
+    double randy = x * f;
+    x ^= x << 13;
+    x ^= x >> 17;
+    x ^= x << 5;
+    double randz = x * f;
+    y = x;
+
+    // calculate coordinates of random point and check its position
+    double px = xmin + (xmax - xmin) * randx;
+    double py = ymin + (ymax - ymin) * randy;
+    double pz = zmin + (zmax - zmin) * randz;
+    icount += contains(px, py, pz);
+  }
+  return (xmax - xmin) * (ymax - ymin) * (zmax -zmin) * icount / npoints;
+}
+
 const GeoShapeUnion & GeoShape::add (const GeoShape& shape) const
 {
   GeoShapeUnion *unionNode = new GeoShapeUnion (this, &shape);
@@ -39,3 +89,138 @@ const GeoShapeShift & GeoShape::operator << (const GeoTrf::Transform3D &shift) c
   GeoShapeShift *shiftNode = new GeoShapeShift (this, shift);
   return *shiftNode;
 }
+
+void GeoShape::diskExtent(double rmin, double rmax, double sphi, double dphi,
+                          double& xmin, double& ymin, double& xmax, double& ymax)
+{
+#ifndef M_PI
+  constexpr double M_PI = 3.14159265358979323846;
+#endif
+  xmin = ymin =-rmax;
+  xmax = ymax = rmax;
+  if (dphi >= 2. * M_PI) return;
+
+  double sinStart = std::sin(sphi);
+  double cosStart = std::cos(sphi);
+  double sinEnd = std::sin(sphi + dphi);
+  double cosEnd = std::cos(sphi + dphi);
+
+  // get start and end quadrants, sixteen possible cases
+  //
+  //      1 | 0
+  //     ---+---
+  //      3 | 2
+  //
+  int icase = (cosEnd < 0) ? 1 : 0;
+  if (sinEnd   < 0) icase += 2;
+  if (cosStart < 0) icase += 4;
+  if (sinStart < 0) icase += 8;
+
+  switch (icase)
+  {
+  // start quadrant 0
+  case  0:                       // start->end : 0->0
+    if (sinEnd < sinStart) break;
+    xmin = rmin*cosEnd;
+    ymin = rmin*sinStart;
+    xmax = rmax*cosStart;
+    ymax = rmax*sinEnd;
+    break;
+  case  1:                       // start->end : 0->1
+    xmin = rmax*cosEnd;
+    ymin = std::min(rmin*sinStart, rmin*sinEnd);
+    xmax = rmax*cosStart;
+    ymax = rmax;
+    break;
+  case  2:                       // start->end : 0->2
+    xmin =-rmax;
+    ymin =-rmax;
+    xmax = std::max(rmax*cosStart, rmax*cosEnd);
+    ymax = rmax;
+    break;
+  case  3:                       // start->end : 0->3
+    xmin =-rmax;
+    ymin = rmax*sinEnd;
+    xmax = rmax*cosStart;
+    ymax = rmax;
+    break;
+  // start quadrant 1
+  case  4:                       // start->end : 1->0
+    xmin =-rmax;
+    ymin =-rmax;
+    xmax = rmax;
+    ymax = std::max(rmax*sinStart, rmax*sinEnd);
+    break;
+  case  5:                       // start->end : 1->1
+    if (sinEnd > sinStart) break;
+    xmin = rmax*cosEnd;
+    ymin = rmin*sinEnd;
+    xmax = rmin*cosStart;
+    ymax = rmax*sinStart;
+    break;
+  case  6:                       // start->end : 1->2
+    xmin =-rmax;
+    ymin =-rmax;
+    xmax = rmax*cosEnd;
+    ymax = rmax*sinStart;
+    break;
+  case  7:                       // start->end : 1->3
+    xmin =-rmax;
+    ymin = rmax*sinEnd;
+    xmax = std::max(rmin*cosStart, rmin*cosEnd);
+    ymax = rmax*sinStart;
+    break;
+  // start quadrant 2
+  case  8:                       // start->end : 2->0
+    xmin = std::min(rmin*cosStart, rmin*cosEnd);
+    ymin = rmax*sinStart;
+    xmax = rmax;
+    ymax = rmax*sinEnd;
+    break;
+  case  9:                       // start->end : 2->1
+    xmin = rmax*cosEnd;
+    ymin = rmax*sinStart;
+    xmax = rmax;
+    ymax = rmax;
+    break;
+  case 10:                       // start->end : 2->2
+    if (sinEnd < sinStart) break;
+    xmin = rmin*cosStart;
+    ymin = rmax*sinStart;
+    xmax = rmax*cosEnd;
+    ymax = rmin*sinEnd;
+    break;
+  case 11:                       // start->end : 2->3
+    xmin =-rmax;
+    ymin = std::min(rmax*sinStart, rmax*sinEnd);
+    xmax = rmax;
+    ymax = rmax;
+    break;
+  // start quadrant 3
+  case 12:                       // start->end : 3->0
+    xmin = rmax*cosStart;
+    ymin =-rmax;
+    xmax = rmax;
+    ymax = rmax*sinEnd;
+    break;
+  case 13:                       // start->end : 3->1
+    xmin = std::min(rmax*cosStart, rmax*cosEnd);
+    ymin =-rmax;
+    xmax = rmax;
+    ymax = rmax;
+    break;
+  case 14:                       // start->end : 3->2
+    xmin = rmax*cosStart;
+    ymin =-rmax;
+    xmax = rmax*cosEnd;
+    ymax = std::max(rmin*sinStart, rmin*sinEnd);
+    break;
+  case 15:                       // start->end : 3->3
+    if (sinEnd > sinStart) break;
+    xmin = rmax*cosStart;
+    ymin = rmax*sinEnd;
+    xmax = rmin*cosEnd;
+    ymax = rmin*sinStart;
+    break;
+  }
+}
diff --git a/GeoModelCore/GeoModelKernel/src/GeoShapeIntersection.cxx b/GeoModelCore/GeoModelKernel/src/GeoShapeIntersection.cxx
index fccc9d3c91713f0a55a175253bf4772fc051d97c..d2f2214761a9170a432dc7a419cc9e9846ed8f8d 100755
--- a/GeoModelCore/GeoModelKernel/src/GeoShapeIntersection.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoShapeIntersection.cxx
@@ -27,11 +27,27 @@ GeoShapeIntersection::~GeoShapeIntersection()
 
 double GeoShapeIntersection::volume () const
 {
-  GeoPolyhedrizeAction a;
-  exec(&a);
-  const GeoPolyhedron *poly = a.getPolyhedron();
-  double vol = poly->GetVolume ();
-  return vol;
+  return (fVolume < 0.) ? (fVolume = GeoShape::volume()) : fVolume;
+}
+
+void GeoShapeIntersection::extent (double& xmin, double& ymin, double& zmin,
+                                   double& xmax, double& ymax, double& zmax) const
+{
+  double xminA, yminA, zminA, xmaxA, ymaxA, zmaxA;
+  double xminB, yminB, zminB, xmaxB, ymaxB, zmaxB;
+  getOpA()->extent(xminA, yminA, zminA, xmaxA, ymaxA, zmaxA);
+  getOpB()->extent(xminB, yminB, zminB, xmaxB, ymaxB, zmaxB);
+  xmin = std::max(xminA, xminB);
+  ymin = std::max(yminA, yminB);
+  zmin = std::max(zminA, zminB);
+  xmax = std::min(xmaxA, xmaxB);
+  ymax = std::min(ymaxA, ymaxB);
+  zmax = std::min(zmaxA, zmaxB);
+}
+
+bool GeoShapeIntersection::contains (double x, double y, double z) const
+{
+  return (getOpA()->contains(x, y, z)) ? getOpB()->contains(x, y, z) : false;
 }
 
 const std::string & GeoShapeIntersection::type () const
diff --git a/GeoModelCore/GeoModelKernel/src/GeoShapeShift.cxx b/GeoModelCore/GeoModelKernel/src/GeoShapeShift.cxx
index 3c49d7702c8276d20ef5077a6abfd7903fda2e3f..ec951461ef08e67f9c37310fd163361802799f5f 100755
--- a/GeoModelCore/GeoModelKernel/src/GeoShapeShift.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoShapeShift.cxx
@@ -25,6 +25,44 @@ double GeoShapeShift::volume () const
   return m_op->volume ();
 }
 
+void GeoShapeShift::extent (double& xmin, double& ymin, double& zmin,
+                            double& xmax, double& ymax, double& zmax) const
+{
+  const GeoShape* shape = getOp();
+  const GeoTrf::Transform3D& trans = getX();
+  double x_min, y_min, z_min, x_max, y_max, z_max;
+  shape->extent(x_min, y_min, z_min, x_max, y_max, z_max);
+  GeoTrf::Vector3D vv[8];
+  vv[0] = trans * GeoTrf::Vector3D(x_min, y_min, z_min);
+  vv[1] = trans * GeoTrf::Vector3D(x_max, y_min, z_min);
+  vv[2] = trans * GeoTrf::Vector3D(x_min, y_max, z_min);
+  vv[3] = trans * GeoTrf::Vector3D(x_max, y_max, z_min);
+  vv[4] = trans * GeoTrf::Vector3D(x_min, y_min, z_max);
+  vv[5] = trans * GeoTrf::Vector3D(x_max, y_min, z_max);
+  vv[6] = trans * GeoTrf::Vector3D(x_min, y_max, z_max);
+  vv[7] = trans * GeoTrf::Vector3D(x_max, y_max, z_max);
+  xmin = xmax = vv[0].x();
+  ymin = ymax = vv[0].y();
+  zmin = zmax = vv[0].z();
+  for (int i = 1; i < 8; ++i)
+  {
+    xmin = std::min(xmin, vv[i].x());
+    ymin = std::min(ymin, vv[i].y());
+    zmin = std::min(zmin, vv[i].z());
+    xmax = std::max(xmax, vv[i].x());
+    ymax = std::max(ymax, vv[i].y());
+    zmax = std::max(zmax, vv[i].z());
+  }
+}
+
+bool GeoShapeShift::contains (double x, double y, double z) const
+{
+  const GeoShape* shape = getOp();
+  const GeoTrf::Transform3D& trans = getX();
+  GeoTrf::Vector3D p = trans.inverse() * GeoTrf::Vector3D(x, y, z);
+  return shape->contains(p.x(), p.y(), p.z());
+}
+
 const std::string & GeoShapeShift::type () const
 {
   return s_classType;
@@ -54,20 +92,20 @@ void GeoShapeShift::exec (GeoShapeAction *action) const
       action->getPath ()->pop ();
       return;
     }
-  
+
   if (action->getDepthLimit ().isValid ()
       && action->getPath ()->getLength () > action->getDepthLimit ())
     {
     }
 
-  else 
-	{
+  else
+        {
     m_op->exec(action);
     if (action->shouldTerminate ())
-	{
+        {
       action->getPath ()->pop ();
       return;
-	}
+        }
   }
   action->getPath()->pop();
 }
diff --git a/GeoModelCore/GeoModelKernel/src/GeoShapeSubtraction.cxx b/GeoModelCore/GeoModelKernel/src/GeoShapeSubtraction.cxx
index 4f3c64f74b902caf7d984b410225121afbc77a35..12c36a1d679ba2aa3d22bb95608819e5e889ccd7 100755
--- a/GeoModelCore/GeoModelKernel/src/GeoShapeSubtraction.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoShapeSubtraction.cxx
@@ -27,11 +27,18 @@ GeoShapeSubtraction::~GeoShapeSubtraction()
 
 double GeoShapeSubtraction::volume () const
 {
-  GeoPolyhedrizeAction a;
-  exec(&a);
-  const GeoPolyhedron *poly = a.getPolyhedron();
-  double vol = poly->GetVolume ();
-  return vol;
+  return (fVolume < 0.) ? (fVolume = GeoShape::volume()) : fVolume;
+}
+
+void GeoShapeSubtraction::extent (double& xmin, double& ymin, double& zmin,
+                                  double& xmax, double& ymax, double& zmax) const
+{
+  getOpA()->extent(xmin, ymin, zmin, xmax, ymax, zmax);
+}
+
+bool GeoShapeSubtraction::contains (double x, double y, double z) const
+{
+  return (!getOpA()->contains(x, y, z)) ? false : !getOpB()->contains(x, y, z);
 }
 
 const std::string & GeoShapeSubtraction::type () const
diff --git a/GeoModelCore/GeoModelKernel/src/GeoShapeUnion.cxx b/GeoModelCore/GeoModelKernel/src/GeoShapeUnion.cxx
index c14c4892c4fd6265494414db76b416351b619d7a..932e637cf026f87e4e70640e2f99aac59271be3f 100755
--- a/GeoModelCore/GeoModelKernel/src/GeoShapeUnion.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoShapeUnion.cxx
@@ -27,11 +27,27 @@ GeoShapeUnion::~GeoShapeUnion()
 
 double GeoShapeUnion::volume () const
 {
-  GeoPolyhedrizeAction a;
-  exec (&a);
-  const GeoPolyhedron *poly = a.getPolyhedron ();
-  double vol = poly->GetVolume ();
-  return vol;
+  return (fVolume < 0.0) ? (fVolume = GeoShape::volume()) : fVolume;
+}
+
+void GeoShapeUnion::extent (double& xmin, double& ymin, double& zmin,
+                            double& xmax, double& ymax, double& zmax) const
+{
+  double xminA, yminA, zminA, xmaxA, ymaxA, zmaxA;
+  double xminB, yminB, zminB, xmaxB, ymaxB, zmaxB;
+  getOpA()->extent(xminA, yminA, zminA, xmaxA, ymaxA, zmaxA);
+  getOpB()->extent(xminB, yminB, zminB, xmaxB, ymaxB, zmaxB);
+  xmin = std::min(xminA, xminB);
+  ymin = std::min(yminA, yminB);
+  zmin = std::min(zminA, zminB);
+  xmax = std::max(xmaxA, xmaxB);
+  ymax = std::max(ymaxA, ymaxB);
+  zmax = std::max(zmaxA, zmaxB);
+}
+
+bool GeoShapeUnion::contains (double x, double y, double z) const
+{
+  return (getOpA()->contains(x, y, z)) ? true : getOpB()->contains(x, y, z);
 }
 
 const std::string & GeoShapeUnion::type () const
@@ -86,4 +102,3 @@ void GeoShapeUnion::exec (GeoShapeAction *action) const
   }
   action->getPath()->pop();
 }
-
diff --git a/GeoModelCore/GeoModelKernel/src/GeoSimplePolygonBrep.cxx b/GeoModelCore/GeoModelKernel/src/GeoSimplePolygonBrep.cxx
index 26adfbb47e31f9c32aefe18126c85a2a3a988647..704c748b28082f0d5969a6f251458e46fe4c1aa2 100755
--- a/GeoModelCore/GeoModelKernel/src/GeoSimplePolygonBrep.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoSimplePolygonBrep.cxx
@@ -30,7 +30,43 @@ double GeoSimplePolygonBrep::volume () const
   {
     area += m_xVertices[k - 1] * m_yVertices[k] - m_xVertices[k] * m_yVertices[k - 1];
   }
-  return fabs(area) * m_dZ;
+  return std::abs(area) * m_dZ;
+}
+
+void GeoSimplePolygonBrep::extent (double& xmin, double& ymin, double& zmin,
+                                   double& xmax, double& ymax, double& zmax) const
+{
+  if (!isValid())
+    throw std::runtime_error ("Extent requested for incomplete simple polygon brep");
+  xmin = xmax = m_xVertices[0];
+  ymin = ymax = m_yVertices[0];
+  for (size_t k = 1; k < getNVertices(); ++k)
+  {
+    double x = m_xVertices[k];
+    double y = m_yVertices[k];
+    xmin = std::min(xmin, x);
+    xmax = std::max(xmax, x);
+    ymin = std::min(ymin, y);
+    ymax = std::max(ymax, y);
+  }
+  zmin =-m_dZ;
+  zmax = m_dZ;
+}
+
+bool GeoSimplePolygonBrep::contains (double x, double y, double z) const
+{
+  if (std::abs(z) - m_dZ > 0.0) return false;
+  size_t nv = getNVertices();
+  bool in = false;
+  for (size_t i = 0, k = nv - 1; i < nv; k = i++)
+  {
+    if ((m_yVertices[i] > y) != (m_yVertices[k] > y))
+    {
+      double ctg = (m_xVertices[k] - m_xVertices[i]) / (m_yVertices[k] - m_yVertices[i]);
+      in ^= (x < (y - m_yVertices[i]) * ctg + m_xVertices[i]);
+    }
+  }
+  return in;
 }
 
 const std::string & GeoSimplePolygonBrep::type() const
@@ -53,4 +89,3 @@ void GeoSimplePolygonBrep::exec(GeoShapeAction *action) const
 {
   action->handleSimplePolygonBrep(this);
 }
-
diff --git a/GeoModelCore/GeoModelKernel/src/GeoTessellatedSolid.cxx b/GeoModelCore/GeoModelKernel/src/GeoTessellatedSolid.cxx
index 44e71c03b6ac3844059af0ec3d566af7ddb8bea6..2d414ae66c2a1a9d35aeff8fca30ae31dbe8bd6a 100755
--- a/GeoModelCore/GeoModelKernel/src/GeoTessellatedSolid.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoTessellatedSolid.cxx
@@ -4,6 +4,7 @@
 
 #include "GeoModelKernel/GeoTessellatedSolid.h"
 #include "GeoModelKernel/GeoShapeAction.h"
+#include <stdexcept>
 
 const std::string GeoTessellatedSolid::s_classType = "TessellatedSolid";
 const ShapeType GeoTessellatedSolid::s_classTypeID = 0x21;
@@ -23,7 +24,7 @@ double GeoTessellatedSolid::volume() const
   if (!isValid ())
     throw std::runtime_error ("Volume requested for incomplete tessellated solid");
   double v = 0.;
-  for (size_t i = 0; i < m_facets.size(); ++i)
+  for (size_t i = 0; i < getNumberOfFacets(); ++i)
   {
     GeoFacet* facet = getFacet(i);
     GeoTrf::Vector3D e1 = facet->getVertex(2) - facet->getVertex(0);
@@ -37,6 +38,41 @@ double GeoTessellatedSolid::volume() const
   return v*(1./6.);
 }
 
+void GeoTessellatedSolid::extent (double& xmin, double& ymin, double& zmin,
+                                  double& xmax, double& ymax, double& zmax) const
+{
+  if (!isValid ())
+    throw std::runtime_error ("Extent requested for incomplete tessellated solid");
+  GeoFacet* facet = getFacet(0);
+  GeoTrf::Vector3D vertex = facet->getVertex(0);
+  xmin = xmax = vertex.x();
+  ymin = ymax = vertex.y();
+  zmin = zmax = vertex.z();
+  for (size_t i = 0; i < getNumberOfFacets(); ++i)
+  {
+    facet = getFacet(i);
+    for (size_t k = 0; k < facet->getNumberOfVertices(); ++k)
+    {
+      vertex = facet->getVertex(k);
+      double x = vertex.x();
+      double y = vertex.y();
+      double z = vertex.z();
+      xmin = std::min(xmin, x);
+      ymin = std::min(ymin, y);
+      zmin = std::min(zmin, z);
+      xmax = std::max(xmax, x);
+      ymax = std::max(ymax, y);
+      zmax = std::max(zmax, z);
+    }
+  }
+}
+
+bool GeoTessellatedSolid::contains (double x, double y, double z) const
+{
+  throw std::runtime_error ("GeoTessellatedSolid::contains(x,y,z) is not implemented");
+  return false;
+}
+
 const std::string& GeoTessellatedSolid::type() const
 {
   return s_classType;
@@ -57,7 +93,7 @@ void GeoTessellatedSolid::addFacet(GeoFacet* facet)
   facet->ref();
   m_facets.push_back(facet);
 }
-  
+
 GeoFacet* GeoTessellatedSolid::getFacet(size_t index) const
 {
   return (index<m_facets.size() ? m_facets[index] : 0);
diff --git a/GeoModelCore/GeoModelKernel/src/GeoTorus.cxx b/GeoModelCore/GeoModelKernel/src/GeoTorus.cxx
index 1a517e361b7264758d91b976032c48c8c26923a1..a9c7cb28f625d7a100442f8c735bd17b8aa5a0ae 100755
--- a/GeoModelCore/GeoModelKernel/src/GeoTorus.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoTorus.cxx
@@ -18,20 +18,54 @@ GeoTorus::GeoTorus (double Rmin, double Rmax, double Rtor, double SPhi, double D
 {
 }
 
-
 GeoTorus::~GeoTorus()
 {
 }
 
-//## Other Operations (implementation)
 double GeoTorus::volume () const
 {
 #ifndef M_PI
-  double M_PI = acos (-1.0);
-#endif 
+  constexpr double M_PI = 3.14159265358979323846;
+#endif
   return m_dPhi * M_PI * m_rTor *(m_rMax * m_rMax - m_rMin * m_rMin);
 }
 
+void GeoTorus::extent (double& xmin, double& ymin, double& zmin,
+                       double& xmax, double& ymax, double& zmax) const
+{
+  double rmin = m_rTor - m_rMax;
+  double rmax = m_rTor + m_rMax;
+  GeoShape::diskExtent(rmin, rmax, m_sPhi, m_dPhi, xmin, ymin, xmax, ymax);
+  zmin =-m_rMax;
+  zmax = m_rMax;
+}
+
+bool GeoTorus::contains (double x, double y, double z) const
+{
+#ifndef M_PI
+  constexpr double M_PI = 3.14159265358979323846;
+#endif
+  double rr = x * x + y * y;
+  if (std::abs(z) - m_rMax > 0.0 || rr == 0.0) return false;
+  double rmin = m_rTor - m_rMax;
+  double rmax = m_rTor + m_rMax;
+  if (rr > rmax * rmax || rr < rmin * rmin) return false;
+
+  double t = (1.0 - m_rTor / std::sqrt(rr));
+  double xt = x * t;
+  double yt = y * t;
+  rr = xt * xt + yt * yt + z * z;
+  if (rr > m_rMax * m_rMax || rr < m_rMin * m_rMin) return false;
+  if (m_dPhi >= 2.0 * M_PI) return true;
+
+  GeoTrf::Vector2D r(x, y);
+  GeoTrf::Vector2D ns(std::sin(m_sPhi), -std::cos(m_sPhi));
+  GeoTrf::Vector2D ne(-std::sin(m_sPhi + m_dPhi), std::cos(m_sPhi + m_dPhi));
+  double ds = ns.dot(r);
+  double de = ne.dot(r);
+  return (m_dPhi <= M_PI) ? (ds <= 0 && de <= 0) : (ds <= 0 || de <= 0);
+}
+
 const std::string & GeoTorus::type () const
 {
   return s_classType;
diff --git a/GeoModelCore/GeoModelKernel/src/GeoTrap.cxx b/GeoModelCore/GeoModelKernel/src/GeoTrap.cxx
index d746707087019fa02a48272addc5dd9fedadf2fd..6f7493d72fb23e0f62ba1e2518f446d156776a24 100755
--- a/GeoModelCore/GeoModelKernel/src/GeoTrap.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoTrap.cxx
@@ -29,15 +29,86 @@ GeoTrap::~GeoTrap()
 
 double GeoTrap::volume () const
 {
-  double fDz  = m_zHalfLength;
-  double fDy1 = m_dydzn;
-  double fDx1 = m_dxdyndzn;
-  double fDx2 = m_dxdypdzn;
-  double fDy2 = m_dydzp;
-  double fDx3 = m_dxdyndzp;
-  double fDx4 = m_dxdypdzp;
-  return fDz * ((fDx1 + fDx2 + fDx3 + fDx4) * (fDy1 + fDy2) +
-                (fDx4 + fDx3 - fDx2 - fDx1) * (fDy2 - fDy1) * (1./3.));
+  double dz = m_zHalfLength;
+  double dy1 = m_dydzn;
+  double dx1 = m_dxdyndzn;
+  double dx2 = m_dxdypdzn;
+  double dy2 = m_dydzp;
+  double dx3 = m_dxdyndzp;
+  double dx4 = m_dxdypdzp;
+  return dz * ((dx1 + dx2 + dx3 + dx4) * (dy1 + dy2) +
+               (dx4 + dx3 - dx2 - dx1) * (dy2 - dy1) * (1./3.));
+}
+
+void GeoTrap::extent (double& xmin, double& ymin, double& zmin,
+                      double& xmax, double& ymax, double& zmax) const
+{
+  double dz = m_zHalfLength;
+  double dy1 = m_dydzn;
+  double dx1 = m_dxdyndzn;
+  double dx2 = m_dxdypdzn;
+  double dy2 = m_dydzp;
+  double dx3 = m_dxdyndzp;
+  double dx4 = m_dxdypdzp;
+  double alpha1 = m_angleydzn;
+  double alpha2 = m_angleydzp;
+
+  double tanTheta = std::tan(m_theta);
+  double dzTthetaCphi = dz * tanTheta * std::cos(m_phi);
+  double dzTthetaSphi = dz * tanTheta * std::sin(m_phi);
+  double dy1Talpha1 = dy1 * std::tan(alpha1);
+  double dy2Talpha2 = dy2 * std::tan(alpha2);
+
+  double x1 =-dzTthetaCphi - dy1Talpha1 - dx1;
+  double x2 =-dzTthetaCphi + dy1Talpha1 - dx2;
+  double x3 = dzTthetaCphi - dy2Talpha2 - dx3;
+  double x4 = dzTthetaCphi + dy2Talpha2 - dx4;
+  double x5 =-dzTthetaCphi - dy1Talpha1 + dx1;
+  double x6 =-dzTthetaCphi + dy1Talpha1 + dx2;
+  double x7 = dzTthetaCphi - dy2Talpha2 + dx3;
+  double x8 = dzTthetaCphi + dy2Talpha2 + dx4;
+
+  double y1 =-dzTthetaSphi - dy1;
+  double y2 = dzTthetaSphi - dy2;
+  double y3 =-dzTthetaSphi + dy1;
+  double y4 = dzTthetaSphi + dy2;
+
+  xmin = std::min(std::min(std::min(x1, x2), x3), x4);
+  xmax = std::max(std::max(std::max(x5, x6), x7), x8);
+  ymin = std::min(y1, y2);
+  ymax = std::max(y3, y4);
+  zmin =-dz;
+  zmax = dz;
+}
+
+bool GeoTrap::contains (double x, double y, double z) const
+{
+  double z0 = z;
+  if (std::abs(z0) - m_zHalfLength > 0.0) return false;
+
+  double dz = m_zHalfLength;
+  double dy1 = m_dydzn;
+  double dx1 = m_dxdyndzn;
+  double dx2 = m_dxdypdzn;
+  double dy2 = m_dydzp;
+  double dx3 = m_dxdyndzp;
+  double dx4 = m_dxdypdzp;
+  double alpha1 = m_angleydzn;
+  double alpha2 = m_angleydzp;
+
+  double tz = 0.5 * (1.0 + z0 / dz);
+  double tanTheta = std::tan(m_theta);
+  double tanAlpha = std::tan(alpha1 + (alpha2 - alpha1) * tz);
+  double dy = dy1 + (dy2 - dy1) * tz;
+  double y0 = y - z0 * tanTheta * std::sin(m_phi);
+  if (std::abs(y0) - dy > 0.0) return false;
+
+  double ty = 0.5 * (1.0 + y0 / dy);
+  double dxneg = dx1 + (dx3 - dx1) * tz;
+  double dxpos = dx2 + (dx4 - dx2) * tz;
+  double dx = dxneg + (dxpos - dxneg) * ty;
+  double x0 = x - y0 * tanAlpha - z0 * tanTheta * std::cos(m_phi);
+  return (std::abs(x0) - dx <= 0.0);
 }
 
 const std::string & GeoTrap::type () const
diff --git a/GeoModelCore/GeoModelKernel/src/GeoTrd.cxx b/GeoModelCore/GeoModelKernel/src/GeoTrd.cxx
index a51167b9f4a5e7c95a070d3f94b0151058580294..51a700e83680af05c4cd5eae2d03709088b001b6 100755
--- a/GeoModelCore/GeoModelKernel/src/GeoTrd.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoTrd.cxx
@@ -30,6 +30,28 @@ double GeoTrd::volume () const
     ((fDx1 + fDx2) * (fDy1 + fDy2) + (fDx2 - fDx1) * (fDy2 - fDy1) * (1./3.));
 }
 
+void GeoTrd::extent (double& xmin, double& ymin, double& zmin,
+                     double& xmax, double& ymax, double& zmax) const
+{
+  xmax = std::max(m_xHalfLength1, m_xHalfLength2);
+  ymax = std::max(m_yHalfLength1, m_yHalfLength2);
+  zmax = m_zHalfLength;
+  xmin = -xmax;
+  ymin = -ymax;
+  zmin = -zmax;
+}
+
+bool GeoTrd::contains (double x, double y, double z) const
+{
+  if (std::abs(z) - m_zHalfLength > 0.0) return false;
+  double t = 0.5 * (1.0 + z / m_zHalfLength);
+  double dx = m_xHalfLength1 + (m_xHalfLength2 - m_xHalfLength1) * t;
+  double dy = m_yHalfLength1 + (m_yHalfLength2 - m_yHalfLength1) * t;
+  double distx = std::abs(x) - dx;
+  double disty = std::abs(y) - dy;
+  return (std::max(distx, disty) <= 0.0);
+}
+
 const std::string & GeoTrd::type () const
 {
   return s_classType;
@@ -44,4 +66,3 @@ void GeoTrd::exec (GeoShapeAction *action) const
 {
   action->handleTrd(this);
 }
-
diff --git a/GeoModelCore/GeoModelKernel/src/GeoTube.cxx b/GeoModelCore/GeoModelKernel/src/GeoTube.cxx
index cdf3c826e2c353088343e6b1703e98cde97b5b6c..7e2c8318089355a9d2d9c0657536960eb90df004 100755
--- a/GeoModelCore/GeoModelKernel/src/GeoTube.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoTube.cxx
@@ -24,9 +24,27 @@ GeoTube::~GeoTube()
 double GeoTube::volume () const
 {
 #ifndef M_PI
-  double M_PI = acos (-1.0);
+  constexpr double M_PI = 3.14159265358979323846;
 #endif
-  return 2 * M_PI * (m_rMax * m_rMax - m_rMin * m_rMin) * m_zHalfLength;
+  return 2.0 * M_PI * (m_rMax * m_rMax - m_rMin * m_rMin) * m_zHalfLength;
+}
+
+void GeoTube::extent (double& xmin, double& ymin, double& zmin,
+                      double& xmax, double& ymax, double& zmax) const
+{
+  xmin =-m_rMax;
+  ymin =-m_rMax;
+  zmin =-m_zHalfLength;
+  xmax = m_rMax;
+  ymax = m_rMax;
+  zmax = m_zHalfLength;
+}
+
+bool GeoTube::contains (double x, double y, double z) const
+{
+  if (std::abs(z) - m_zHalfLength > 0.0) return false;
+  double rr = x * x + y * y;
+  return ((rr <= m_rMax * m_rMax) && (rr >= m_rMin * m_rMin));
 }
 
 const std::string & GeoTube::type () const
diff --git a/GeoModelCore/GeoModelKernel/src/GeoTubs.cxx b/GeoModelCore/GeoModelKernel/src/GeoTubs.cxx
index 1ffaf1d7dc3c34cb95d6429099a73a8bb2331ecb..8568d5dd9d53a45e8b723d62f9764d53e9c496df 100755
--- a/GeoModelCore/GeoModelKernel/src/GeoTubs.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoTubs.cxx
@@ -4,6 +4,7 @@
 
 #include "GeoModelKernel/GeoTubs.h"
 #include "GeoModelKernel/GeoShapeAction.h"
+#include <cmath>
 
 const std::string GeoTubs::s_classType = "Tubs";
 const ShapeType GeoTubs::s_classTypeID = 0x18;
@@ -26,6 +27,32 @@ double GeoTubs::volume () const
   return m_dPhi * (m_rMax * m_rMax - m_rMin * m_rMin) * m_zHalfLength;
 }
 
+void GeoTubs::extent (double& xmin, double& ymin, double& zmin,
+                      double& xmax, double& ymax, double& zmax) const
+{
+  GeoShape::diskExtent(m_rMin, m_rMax, m_sPhi, m_dPhi, xmin, ymin, xmax, ymax);
+  zmin =-m_zHalfLength;
+  zmax = m_zHalfLength;
+}
+
+bool GeoTubs::contains (double x, double y, double z) const
+{
+#ifndef M_PI
+  constexpr double M_PI = 3.14159265358979323846;
+#endif
+  if (std::abs(z) - m_zHalfLength > 0.0) return false;
+  double rr = x * x + y * y;
+  if (rr > m_rMax * m_rMax || rr < m_rMin * m_rMin) return false;
+  if (m_dPhi >= 2.0 * M_PI) return true;
+
+  GeoTrf::Vector2D r(x, y);
+  GeoTrf::Vector2D ns(std::sin(m_sPhi), -std::cos(m_sPhi));
+  GeoTrf::Vector2D ne(-std::sin(m_sPhi + m_dPhi), std::cos(m_sPhi + m_dPhi));
+  double ds = ns.dot(r);
+  double de = ne.dot(r);
+  return (m_dPhi <= M_PI) ? (ds <= 0 && de <= 0) : (ds <= 0 || de <= 0);
+}
+
 const std::string & GeoTubs::type () const
 {
   return s_classType;
@@ -40,4 +67,3 @@ void GeoTubs::exec (GeoShapeAction *action) const
 {
   action->handleTubs(this);
 }
-
diff --git a/GeoModelCore/GeoModelKernel/src/GeoTwistedTrap.cxx b/GeoModelCore/GeoModelKernel/src/GeoTwistedTrap.cxx
index c96bbaa6661a7d5535be1f7e2c0372255e324a01..3054709214c9f2f35fe89a6c59b23c23f5041dc1 100755
--- a/GeoModelCore/GeoModelKernel/src/GeoTwistedTrap.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoTwistedTrap.cxx
@@ -4,11 +4,12 @@
 
 #include "GeoModelKernel/GeoTwistedTrap.h"
 #include "GeoModelKernel/GeoShapeAction.h"
+#include <cmath>
 #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,
@@ -16,16 +17,16 @@ GeoTwistedTrap::GeoTwistedTrap(double  pPhiTwist,
                                double  pDy,
                                double  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_dz(pDz),
-m_alph(0.),
-m_phiTwist(pPhiTwist)
+  m_phi(0.),
+  m_dy1(pDy),
+  m_dx1(pDx1),
+  m_dx2(pDx2),
+  m_dy2(pDy),
+  m_dx3(pDx1),
+  m_dx4(pDx2),
+  m_dz(pDz),
+  m_alph(0.),
+  m_phiTwist(pPhiTwist)
 {
 }
 
@@ -42,16 +43,16 @@ GeoTwistedTrap(double  pPhiTwist,  // twist angle
                double  pDx4,    // half x length at +pDz,+pDy
                double  pAlph )  // tilt angle
 : 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_dz(pDz),
-m_alph(pAlph),
-m_phiTwist(pPhiTwist)
+  m_phi(pPhi),
+  m_dy1(pDy1),
+  m_dx1(pDx1),
+  m_dx2(pDx2),
+  m_dy2(pDy2),
+  m_dx3(pDx3),
+  m_dx4(pDx4),
+  m_dz(pDz),
+  m_alph(pAlph),
+  m_phiTwist(pPhiTwist)
 {
      if  ( ! ( ( m_dx1  > 2*m_CarTolerance)
             && ( m_dx2  > 2*m_CarTolerance)
@@ -60,9 +61,9 @@ m_phiTwist(pPhiTwist)
             && ( 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 )
+            && ( std::abs(m_phiTwist) > 2*m_AngTolerance )
+            && ( std::abs(m_phiTwist) < SYSTEM_OF_UNITS::pi/2 )
+            && ( std::abs(m_alph) < SYSTEM_OF_UNITS::pi/2 )
             && ( m_theta < SYSTEM_OF_UNITS::pi/2 && m_theta >= 0 ) )
          )
      {
@@ -79,7 +80,7 @@ m_phiTwist(pPhiTwist)
                << " 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;
-    
+
      }
 }
 
@@ -94,6 +95,65 @@ double GeoTwistedTrap::volume () const
                  (m_dx4 + m_dx3 - m_dx2 - m_dx1) * (m_dy2 - m_dy1) * (1./3.));
 }
 
+void GeoTwistedTrap::extent (double& xmin, double& ymin, double& zmin,
+                             double& xmax, double& ymax, double& zmax) const
+{
+  double cosPhi = std::cos(m_phi);
+  double sinPhi = std::sin(m_phi);
+  double tanTheta = std::tan(m_theta);
+  double tanAlpha = std::tan(m_alph);
+
+  double xmid1 = m_dy1 * tanAlpha;
+  double x1 = std::abs(xmid1 + m_dx1);
+  double x2 = std::abs(xmid1 - m_dx1);
+  double x3 = std::abs(xmid1 + m_dx2);
+  double x4 = std::abs(xmid1 - m_dx2);
+  double xmax1 = std::max(std::max(std::max(x1, x2), x3), x4);
+  double rmax1 = std::hypot(xmax1, m_dy1);
+
+  double xmid2 = m_dy2 * tanAlpha;
+  double x5 = std::abs(xmid2 + m_dx3);
+  double x6 = std::abs(xmid2 - m_dx3);
+  double x7 = std::abs(xmid2 + m_dx4);
+  double x8 = std::abs(xmid2 - m_dx4);
+  double xmax2 = std::max(std::max(std::max(x5, x6), x7), x8);
+  double rmax2 = hypot(xmax2, m_dy2);
+
+  double x0 = m_dz * tanTheta * cosPhi;
+  double y0 = m_dz * tanTheta * sinPhi;
+  xmin = std::min(-x0 - rmax1, x0 - rmax2);
+  ymin = std::min(-y0 - rmax1, y0 - rmax2);
+  xmax = std::max(-x0 + rmax1, x0 + rmax2);
+  ymax = std::max(-y0 + rmax1, y0 + rmax2);
+  zmin =-m_dz;
+  zmax = m_dz;
+}
+
+bool GeoTwistedTrap::contains (double x, double y, double z) const
+{
+  double z0 = z;
+  if (std::abs(z0) - m_dz > 0.0) return false;
+  double tz = 0.5 * (1.0 + z0 / m_dz);
+
+  double twist = -0.5 * m_phiTwist + m_phiTwist * tz;
+  double cosTwist = std::cos(-twist);
+  double sinTwist = std::sin(-twist);
+  double tanTheta = std::tan(m_theta);
+  double xc = z0 * tanTheta * std::cos(m_phi);
+  double yc = z0 * tanTheta * std::sin(m_phi);
+
+  double dy = m_dy1 + (m_dy2 - m_dy1) * tz;
+  double y0 = sinTwist * (x - xc) + cosTwist * (y - yc);
+  if (std::abs(y0) - dy > 0.0) return false;
+  double ty = 0.5 * (1.0 + y0 / dy);
+
+  double dxneg = m_dx1 + (m_dx3 - m_dx1) * tz;
+  double dxpos = m_dx2 + (m_dx4 - m_dx2) * tz;
+  double dx = dxneg + (dxpos - dxneg) * ty;
+  double x0 = cosTwist * (x - xc) - sinTwist * (y - yc) - y0 * std::tan(m_alph);
+  return (std::abs(x0) - dx <= 0.0);
+}
+
 const std::string & GeoTwistedTrap::type () const
 {
   return s_classType;
@@ -104,7 +164,7 @@ ShapeType GeoTwistedTrap::typeID () const
   return s_classTypeID;
 }
 
-void GeoTwistedTrap::exec (GeoShapeAction *action) const
+void GeoTwistedTrap::exec (GeoShapeAction* action) const
 {
-  action->handleTwistedTrap(this); 
+  action->handleTwistedTrap(this);
 }
diff --git a/GeoModelCore/GeoModelKernel/src/GeoUnidentifiedShape.cxx b/GeoModelCore/GeoModelKernel/src/GeoUnidentifiedShape.cxx
index cdae19b1b96f18adba5d4a98350b4e2ac533be8e..a5292a185b981b0fdc4fec9e730c98e5ed30c02e 100644
--- a/GeoModelCore/GeoModelKernel/src/GeoUnidentifiedShape.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoUnidentifiedShape.cxx
@@ -1,6 +1,7 @@
 #include "GeoModelKernel/GeoUnidentifiedShape.h"
 #include "GeoModelKernel/Query.h"
 #include "GeoModelKernel/GeoShapeAction.h"
+#include <stdexcept>
 
 const std::string         GeoUnidentifiedShape::_classType="UnidentifiedShape";
 const ShapeType           GeoUnidentifiedShape::_classTypeID=0xFFFFFFFF;
@@ -18,7 +19,7 @@ GeoUnidentifiedShape::GeoUnidentifiedShape(const std::string & name, const std::
   _name(name),
   _asciiData(asciiData) {}
 
-// Constructor with volume: 
+// Constructor with volume:
 GeoUnidentifiedShape::GeoUnidentifiedShape(const std::string & name, const std::string & asciiData, double volume):
   _name(name),
   _asciiData(asciiData),
@@ -38,7 +39,22 @@ const std::string & GeoUnidentifiedShape::asciiData() const {
 double GeoUnidentifiedShape::volume () const {
   return _volume;
 }
-  
+
+// Returns the bonding box of the shape
+void GeoUnidentifiedShape::extent (double& xmin, double& ymin, double& zmin,
+                                   double& xmax, double& ymax, double& zmax) const
+{
+  throw std::runtime_error ("GeoUndefinedShape::extent is not implemented");
+  xmin = ymin = zmin = xmax = ymax = zmax = 0.;
+}
+
+// Returns true if the shape contains the point, false otherwise
+bool GeoUnidentifiedShape::contains (double x, double y, double z) const
+{
+  throw std::runtime_error ("GeoUndefinedShape::contains(x,y,z) is not implemented");
+  return false;
+}
+
 // Returns the shape type, as a string.
 const std::string & GeoUnidentifiedShape::type () const {
   return _classType;