From ad3beeb7d8e05e411353b672a1e3ac50a4b6d2cd Mon Sep 17 00:00:00 2001
From: Andrei Gheata <Andrei.Gheata@cern.ch>
Date: Mon, 10 Dec 2018 17:29:38 +0100
Subject: [PATCH] Added GDML support for: opticalsurface, skinsurface,
 bordersurface, matrix, property, including test macro

---
 geom/gdml/inc/TGDMLParse.h               |  19 +-
 geom/gdml/inc/TGDMLWrite.h               |  14 +
 geom/gdml/src/TGDMLParse.cxx             | 707 ++++++++++++++++-------
 geom/gdml/src/TGDMLWrite.cxx             | 162 +++++-
 geom/geom/CMakeLists.txt                 |   4 +-
 geom/geom/inc/LinkDef1.h                 |   7 +
 geom/geom/inc/TGDMLMatrix.h              |  61 ++
 geom/geom/inc/TGeoManager.h              |  26 +-
 geom/geom/inc/TGeoMaterial.h             |  13 +-
 geom/geom/inc/TGeoOpticalSurface.h       | 207 +++++++
 geom/geom/src/TGDMLMatrix.cxx            |  99 ++++
 geom/geom/src/TGeoManager.cxx            |  99 ++++
 geom/geom/src/TGeoMaterial.cxx           |  29 +
 geom/geom/src/TGeoOpticalSurface.cxx     | 300 ++++++++++
 tutorials/geom/gdml/opticalsurfaces.gdml |  78 +++
 tutorials/geom/gdml/testoptical.C        |  69 +++
 16 files changed, 1679 insertions(+), 215 deletions(-)
 create mode 100644 geom/geom/inc/TGDMLMatrix.h
 create mode 100644 geom/geom/inc/TGeoOpticalSurface.h
 create mode 100644 geom/geom/src/TGDMLMatrix.cxx
 create mode 100644 geom/geom/src/TGeoOpticalSurface.cxx
 create mode 100644 tutorials/geom/gdml/opticalsurfaces.gdml
 create mode 100644 tutorials/geom/gdml/testoptical.C

diff --git a/geom/gdml/inc/TGDMLParse.h b/geom/gdml/inc/TGDMLParse.h
index 9bba43c8d31..5f198564d72 100644
--- a/geom/gdml/inc/TGDMLParse.h
+++ b/geom/gdml/inc/TGDMLParse.h
@@ -24,6 +24,8 @@
 #include <vector>
 #include <iostream>
 
+class TGDMLMatrix;
+
 /*************************************************************************
  * TGDMLRefl - helper class for the import of GDML to ROOT.              *
  *************************************************************************/
@@ -137,6 +139,7 @@ private:
    double            Evaluate(const char* evalline);
    const char*       NameShort(const char* name);
    double            Value(const char *svalue) const;
+   void              DefineConstants();
 
    //'define' section
    XMLNodePointer_t  ConProcess(TXMLEngine* gdml, XMLNodePointer_t node, XMLAttrPointer_t attr);
@@ -144,15 +147,15 @@ private:
    XMLNodePointer_t  QuantityProcess(TXMLEngine* gdml, XMLNodePointer_t node, XMLAttrPointer_t attr);
    XMLNodePointer_t  RotProcess(TXMLEngine* gdml, XMLNodePointer_t node, XMLAttrPointer_t attr);
    XMLNodePointer_t  SclProcess(TXMLEngine* gdml, XMLNodePointer_t node, XMLAttrPointer_t attr);
+   XMLNodePointer_t  MatrixProcess(TXMLEngine* gdml, XMLNodePointer_t node, XMLAttrPointer_t attr);
 
    //'materials' section
-  XMLNodePointer_t  IsoProcess(TXMLEngine* gdml, XMLNodePointer_t node, XMLNodePointer_t parentn);
-  XMLNodePointer_t  EleProcess(TXMLEngine* gdml, XMLNodePointer_t node, XMLNodePointer_t parentn, Bool_t hasIsotopes, Bool_t hasIsotopesExtended);
-  //XMLNodePointer_t  EleProcess(TXMLEngine* gdml, XMLNodePointer_t node, XMLNodePointer_t parentn);
-  XMLNodePointer_t  MatProcess(TXMLEngine* gdml, XMLNodePointer_t node, XMLAttrPointer_t attr,  int z);
-  //XMLNodePointer_t  MatProcess(TXMLEngine* gdml, XMLNodePointer_t node, XMLAttrPointer_t attr);
+   XMLNodePointer_t  IsoProcess(TXMLEngine* gdml, XMLNodePointer_t node, XMLNodePointer_t parentn);
+   XMLNodePointer_t  EleProcess(TXMLEngine* gdml, XMLNodePointer_t node, XMLNodePointer_t parentn, Bool_t hasIsotopes, Bool_t hasIsotopesExtended);
+   XMLNodePointer_t  MatProcess(TXMLEngine* gdml, XMLNodePointer_t node, XMLAttrPointer_t attr,  int z);
 
    //'solids' section
+   XMLNodePointer_t  OpticalSurfaceProcess(TXMLEngine* gdml, XMLNodePointer_t node, XMLAttrPointer_t attr);
    XMLNodePointer_t  BooSolid(TXMLEngine* gdml, XMLNodePointer_t node, XMLAttrPointer_t attr, int num);
    XMLNodePointer_t  Box(TXMLEngine* gdml, XMLNodePointer_t node, XMLAttrPointer_t attr);
    XMLNodePointer_t  Paraboloid(TXMLEngine* gdml, XMLNodePointer_t node, XMLAttrPointer_t attr);
@@ -180,6 +183,8 @@ private:
    XMLNodePointer_t  VolProcess(TXMLEngine* gdml, XMLNodePointer_t node);
    XMLNodePointer_t  AssProcess(TXMLEngine* gdml, XMLNodePointer_t node);
    XMLNodePointer_t  UsrProcess(TXMLEngine* gdml, XMLNodePointer_t node);
+   XMLNodePointer_t  SkinSurfaceProcess(TXMLEngine* gdml, XMLNodePointer_t node, XMLAttrPointer_t attr);
+   XMLNodePointer_t  BorderSurfaceProcess(TXMLEngine* gdml, XMLNodePointer_t node, XMLAttrPointer_t attr);
    Int_t             SetAxis(const char* axisString); //Set Axis for Division
 
    //'setup' section
@@ -196,6 +201,8 @@ private:
 
    typedef TGDMMapHelper<TGeoShape> SolMap;
    typedef TGDMMapHelper<TGeoVolume> VolMap;
+   typedef TGDMMapHelper<TGeoNode> PvolMap;
+   typedef TGDMMapHelper<TGDMLMatrix> MatrixMap;
    typedef TGDMMapHelper<TGDMLRefl> ReflSolidMap;
    typedef TGDMMapHelper<const char> FileMap;
    typedef std::map<std::string, std::string> ReflectionsMap;
@@ -213,11 +220,13 @@ private:
    MixMap fmixmap;                //!Map containing mixture names and the TGeoMixture for it
    SolMap fsolmap;                //!Map containing solid names and the TGeoShape for it
    VolMap fvolmap;                //!Map containing volume names and the TGeoVolume for it
+   PvolMap fpvolmap;              //!Map containing placed volume names and the TGeoNode for it
    ReflectionsMap freflectmap;    //!Map containing reflection names and the Solid name ir references to
    ReflSolidMap freflsolidmap;    //!Map containing reflection names and the TGDMLRefl for it - containing refl matrix
    ReflVolMap freflvolmap;        //!Map containing reflected volume names and the solid ref for it
    FileMap ffilemap;              //!Map containing files parsed during entire parsing, with their world volume name
    ConstMap fconsts;              //!Map containing values of constants declared in the file
+   MatrixMap fmatrices;           //!Map containing matrices defined in the GDML file
 
    ClassDef(TGDMLParse, 0)    //imports GDML using DOM and binds it to ROOT
 };
diff --git a/geom/gdml/inc/TGDMLWrite.h b/geom/gdml/inc/TGDMLWrite.h
index 5d16a18c46a..7b19e943562 100644
--- a/geom/gdml/inc/TGDMLWrite.h
+++ b/geom/gdml/inc/TGDMLWrite.h
@@ -34,6 +34,7 @@
 #include "TGeoCompositeShape.h"
 #include "TGeoScaledShape.h"
 #include "TGeoManager.h"
+#include "TGDMLMatrix.h"
 
 #include <map>
 #include <vector>
@@ -47,6 +48,10 @@
 //                                                                        //
 ////////////////////////////////////////////////////////////////////////////
 
+class TGeoOpticalSurface;
+class TGeoSkinSurface;
+class TGeoBorderSurface;
+
 class TGDMLWrite : public TObject {
 public:
    TGDMLWrite();
@@ -134,6 +139,10 @@ private:
    XMLNodePointer_t ExtractMaterials(TList* materialsLst); //result <materials>...
    TString          ExtractSolid(TGeoShape* volShape);     //adds <shape> to <solids>
    void             ExtractVolumes(TGeoVolume* volume);    //result <volume> node...  + corresp. shape
+   void             ExtractMatrices(TObjArray *matrices);  //adds <matrix> to <define>
+   void             ExtractOpticalSurfaces(TObjArray *surfaces); //adds <opticalsurface> to <solids>
+   void             ExtractSkinSurfaces(TObjArray *surfaces);    //adds <skinsurface> to <structure>
+   void             ExtractBorderSurfaces(TObjArray *surfaces);  //adds <bordersurface> to <structure>
 
    // Combined implementation to extract GDML information from the geometry tree
    void WriteGDMLfile(TGeoManager * geomanager, TGeoVolume* volume, TList* materialsLst, const char* filename, TString option);
@@ -142,6 +151,7 @@ private:
    XMLNodePointer_t CreateAtomN(Double_t atom, const char * unit = "g/mole");
    XMLNodePointer_t CreateDN(Double_t density, const char * unit = "g/cm3");
    XMLNodePointer_t CreateFractionN(Double_t percentage, const char * refName);
+   XMLNodePointer_t CreatePropertyN(TNamed const &property);
 
    XMLNodePointer_t CreateIsotopN(TGeoIsotope * isotope, const char * name);
    XMLNodePointer_t CreateElementN(TGeoElement * element, XMLNodePointer_t materials, const char * name);
@@ -175,6 +185,9 @@ private:
    XMLNodePointer_t CreateXtrusionN(TGeoXtru * geoShape);
    XMLNodePointer_t CreateEllipsoidN(TGeoCompositeShape * geoShape, TString elName);
    XMLNodePointer_t CreateElConeN(TGeoScaledShape * geoShape);
+   XMLNodePointer_t CreateOpticalSurfaceN(TGeoOpticalSurface * geoSurf);
+   XMLNodePointer_t CreateSkinSurfaceN(TGeoSkinSurface * geoSurf);
+   XMLNodePointer_t CreateBorderSurfaceN(TGeoBorderSurface * geoSurf);
 
    XMLNodePointer_t CreateCommonBoolN(TGeoCompositeShape *geoShape);
 
@@ -192,6 +205,7 @@ private:
    //nodes to create position, rotation and similar types first-position/rotation...
    XMLNodePointer_t CreatePositionN(const char * name, Xyz position, const char * type = "position", const char * unit = "cm");
    XMLNodePointer_t CreateRotationN(const char * name, Xyz rotation, const char * type = "rotation", const char * unit = "deg");
+   XMLNodePointer_t CreateMatrixN(TGDMLMatrix const *matrix);
    TGeoCompositeShape* CreateFakeCtub(TGeoCtub * geoShape);  //create fake cut tube as intersection
 
    //check name (2nd parameter) whether it is in the list (1st parameter)
diff --git a/geom/gdml/src/TGDMLParse.cxx b/geom/gdml/src/TGDMLParse.cxx
index 1c87a32be34..8c0cf243c31 100644
--- a/geom/gdml/src/TGDMLParse.cxx
+++ b/geom/gdml/src/TGDMLParse.cxx
@@ -77,6 +77,9 @@ When most solids or volumes are added to the geometry they
 
 */
 
+#include "TGDMLParse.h"
+#include "TGDMLMatrix.h"
+
 #include "TGeoManager.h"
 #include "TGeoMatrix.h"
 #include "TXMLEngine.h"
@@ -109,9 +112,12 @@ When most solids or volumes are added to the geometry they
 #include "TGeoShape.h"
 #include "TGeoCompositeShape.h"
 #include "TGeoRegion.h"
-#include "TGDMLParse.h"
+#include "TGeoOpticalSurface.h"
+#include "TGeoSystemOfUnits.h"
+
 #include <stdlib.h>
 #include <string>
+#include <sstream>
 #include <locale>
 
 ClassImp(TGDMLParse);
@@ -160,6 +166,7 @@ TGeoVolume* TGDMLParse::GDMLReadFile(const char* filename)
 
 const char* TGDMLParse::ParseGDML(TXMLEngine* gdml, XMLNodePointer_t node)
 {
+   DefineConstants();
    XMLAttrPointer_t attr = gdml->GetFirstAttr(node);
    const char* name = gdml->GetNodeName(node);
    XMLNodePointer_t parentn = gdml->GetParent(node);
@@ -171,6 +178,7 @@ const char* TGDMLParse::ParseGDML(TXMLEngine* gdml, XMLNodePointer_t node)
    const char* consstr = "constant";
    const char* varistr = "variable";
    const char* quanstr = "quantity";
+   const char* matrstr = "matrix";
    const char* rotastr = "rotation";
    const char* scalstr = "scale";
    const char* elemstr = "element";
@@ -178,7 +186,7 @@ const char* TGDMLParse::ParseGDML(TXMLEngine* gdml, XMLNodePointer_t node)
    const char* matestr = "material";
    const char* volustr = "volume";
    const char* assestr = "assembly";
-   const char* twtrstr = "twistedtrap"; //name changed according to schema
+   const char* twtrstr = "twistedtrap";
    const char* cutTstr = "cutTube";
    const char* bboxstr = "box";
    const char* xtrustr = "xtru";
@@ -202,6 +210,9 @@ const char* TGDMLParse::ParseGDML(TXMLEngine* gdml, XMLNodePointer_t node)
    const char* reflstr = "reflectedSolid";
    const char* ellistr = "ellipsoid";
    const char* elcnstr = "elcone";
+   const char* optsstr = "opticalsurface";
+   const char* skinstr = "skinsurface";
+   const char* bordstr = "bordersurface";
    const char* usrstr  = "userinfo";
    Bool_t hasIsotopes;
    Bool_t hasIsotopesExtended;
@@ -220,6 +231,14 @@ const char* TGDMLParse::ParseGDML(TXMLEngine* gdml, XMLNodePointer_t node)
       node = ConProcess(gdml, node, attr);
    } else if ((strcmp(name, quanstr)) == 0) {
       node = QuantityProcess(gdml, node, attr);
+   } else if ((strcmp(name, matrstr)) == 0) {
+      node = MatrixProcess(gdml, node, attr);
+   } else if ((strcmp(name, optsstr)) == 0) {
+      node = OpticalSurfaceProcess(gdml, node, attr);
+   } else if ((strcmp(name, skinstr)) == 0) {
+      node = SkinSurfaceProcess(gdml, node, attr);
+   } else if ((strcmp(name, bordstr)) == 0) {
+      node = BorderSurfaceProcess(gdml, node, attr);
    }
    //*************eleprocess********************************
 
@@ -429,6 +448,44 @@ XMLNodePointer_t TGDMLParse::ConProcess(TXMLEngine* gdml, XMLNodePointer_t node,
    return node;
 }
 
+
+////////////////////////////////////////////////////////////////////////////////
+/// Define constant expressions used.
+void TGDMLParse::DefineConstants()
+{
+   // Units used in TGeo. Note that they are based on cm/degree/GeV and they are different from Geant4
+   fconsts["mm"] = TGeoUnit::mm;
+   fconsts["millimeter"] = TGeoUnit::mm;
+   fconsts["cm"] = TGeoUnit::cm;
+   fconsts["centimeter"] = TGeoUnit::cm;
+   fconsts["m"] = TGeoUnit::m;
+   fconsts["meter"] = TGeoUnit::m;
+   fconsts["km"] = TGeoUnit::km;
+   fconsts["kilometer"] = TGeoUnit::km;
+   fconsts["rad"] = TGeoUnit::rad;
+   fconsts["radian"] = TGeoUnit::rad;
+   fconsts["deg"] = TGeoUnit::deg;
+   fconsts["degree"] = TGeoUnit::deg;
+   fconsts["pi"] = TGeoUnit::pi;
+   fconsts["twopi"] = TGeoUnit::twopi;
+   fconsts["avogadro"] = TMath::Na();
+   fconsts["gev"] = TGeoUnit::GeV;
+   fconsts["GeV"] = TGeoUnit::GeV;
+   fconsts["mev"] = TGeoUnit::MeV;
+   fconsts["MeV"] = TGeoUnit::MeV;
+   fconsts["kev"] = TGeoUnit::keV;
+   fconsts["keV"] = TGeoUnit::keV;
+   fconsts["ev"] = TGeoUnit::eV;
+   fconsts["eV"] = TGeoUnit::eV;
+   fconsts["s"] = TGeoUnit::s;
+   fconsts["ms"] = TGeoUnit::ms;
+   fconsts["ns"] = TGeoUnit::ns;
+   fconsts["us"] = TGeoUnit::us;
+   fconsts["kg"] = TGeoUnit::kg;
+   fconsts["g"] = TGeoUnit::g;
+   fconsts["mg"] = TGeoUnit::mg;
+}
+
 ////////////////////////////////////////////////////////////////////////////////
 /// In the define section of the GDML file, quantities can be declared.
 /// These are treated the same as constants, but the unit has to be multiplied
@@ -461,6 +518,115 @@ XMLNodePointer_t TGDMLParse::QuantityProcess(TXMLEngine* gdml, XMLNodePointer_t
    return node;
 }
 
+////////////////////////////////////////////////////////////////////////////////
+/// In the define section of the GDML file, matrices
+/// These are referenced by other GDML tags, such as optical surfaces
+XMLNodePointer_t TGDMLParse::MatrixProcess(TXMLEngine* gdml, XMLNodePointer_t node, XMLAttrPointer_t attr)
+{
+   TString name = "";
+   Int_t coldim  = 0;
+   std::string values;
+   TString tempattr;
+
+   while (attr != 0) {
+      tempattr = gdml->GetAttrName(attr);
+      tempattr.ToLower();
+
+      if (tempattr == "name") {
+         name = gdml->GetAttrValue(attr);
+      }
+      if (tempattr == "coldim") {
+         coldim = (Int_t)Value(gdml->GetAttrValue(attr));
+      }
+      if (tempattr == "values") {
+         values = gdml->GetAttrValue(attr);
+      }
+      attr = gdml->GetNextAttr(attr);
+   }
+
+   // Parse the values and create the matrix
+   std::stringstream valueStream(values);
+   std::vector<Double_t> valueList;
+   while (!valueStream.eof())
+   {
+      std::string matrixValue;
+      valueStream >> matrixValue;
+
+      valueList.push_back(Value(matrixValue.c_str()));
+   }
+
+   TGDMLMatrix *matrix = new TGDMLMatrix(name, valueList.size()/coldim, coldim);
+   matrix->SetMatrixAsString(values.c_str());
+   for (size_t i=0; i<valueList.size(); ++i)
+      matrix->Set(i/coldim, i%coldim, valueList[i]);
+
+   gGeoManager->AddGDMLMatrix(matrix);
+   fmatrices[name.Data()] = matrix;
+
+   return node;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+/// In the solids section of the GDML file, optical surfaces can be defined
+///
+XMLNodePointer_t TGDMLParse::OpticalSurfaceProcess(TXMLEngine* gdml, XMLNodePointer_t node, XMLAttrPointer_t attr)
+{
+   TString name, propname, ref;
+   TGeoOpticalSurface::ESurfaceModel model = TGeoOpticalSurface::kMglisur;
+   TGeoOpticalSurface::ESurfaceFinish finish = TGeoOpticalSurface::kFpolished;
+   TGeoOpticalSurface::ESurfaceType type = TGeoOpticalSurface::kTdielectric_metal;
+   Double_t value = 0;
+   TString tempattr;
+
+   while (attr != 0) {
+      tempattr = gdml->GetAttrName(attr);
+      tempattr.ToLower();
+
+      if (tempattr == "name") {
+         name = gdml->GetAttrValue(attr);
+      }
+      if (tempattr == "model") {
+         model = TGeoOpticalSurface::StringToModel(gdml->GetAttrValue(attr));
+      }
+      if (tempattr == "finish") {
+         finish = TGeoOpticalSurface::StringToFinish(gdml->GetAttrValue(attr));
+      }
+      if (tempattr == "type") {
+         type = TGeoOpticalSurface::StringToType(gdml->GetAttrValue(attr));
+      }
+      if (tempattr == "value") {
+         value = Value(gdml->GetAttrValue(attr));
+      }
+      attr = gdml->GetNextAttr(attr);
+   }
+
+   TGeoOpticalSurface *surf = new TGeoOpticalSurface(name, model, finish, type, value);
+
+   XMLNodePointer_t child = gdml->GetChild(node);
+   while (child != 0) {
+      attr = gdml->GetFirstAttr(child);
+      if ((strcmp(gdml->GetNodeName(child), "property")) == 0) {
+         while (attr != 0) {
+            tempattr = gdml->GetAttrName(attr);
+            tempattr.ToLower();
+            if (tempattr == "name") {
+               propname = gdml->GetAttrValue(attr);
+            } else if (tempattr == "ref") {
+               ref = gdml->GetAttrValue(attr);
+               TGDMLMatrix *matrix = fmatrices[ref.Data()];
+               if (!matrix)
+                  Error("OpticalSurfaceProcess", "Reference matrix %s for optical surface %s not found", ref.Data(), name.Data());
+               surf->AddProperty(propname, ref);
+            }
+            attr = gdml->GetNextAttr(attr);
+         }
+      } // loop on child attributes
+      child = gdml->GetNext(child);
+   } // loop on children
+   gGeoManager->AddOpticalSurface(surf);
+   return child;
+}
+
 ////////////////////////////////////////////////////////////////////////////////
 /// Throughout the GDML file, a unit can de specified.   Whether it be
 /// angular or linear, values can be used as well as abbreviations such as
@@ -1074,243 +1240,381 @@ XMLNodePointer_t TGDMLParse::EleProcess(TXMLEngine* gdml, XMLNodePointer_t node,
 
 XMLNodePointer_t TGDMLParse::MatProcess(TXMLEngine* gdml, XMLNodePointer_t node,   XMLAttrPointer_t attr,  int z)
 {
-  //!Map to hold fractions while being processed
-  typedef FracMap::iterator fractions;
+   //!Map to hold fractions while being processed
+   typedef FracMap::iterator fractions;
 //  typedef FracMap::iterator i;
-  FracMap fracmap;
+   FracMap fracmap;
+
+   TGeoManager* mgr = gGeoManager;
+   TGeoElementTable* tab_ele = mgr->GetElementTable();
+   TList properties;
+   properties.SetOwner();
+   // We have to assume the media are monotonic increasing starting with 1
+   static int medid = mgr->GetListOfMedia()->GetSize()+1;
+   XMLNodePointer_t child = gdml->GetChild(node);
+   TString tempattr = "";
+   Int_t ncompo = 0, mixflag = 2;
+   Double_t density = 0;
+   TString name = "";
+   TGeoMixture* mix = 0;
+   TGeoMaterial* mat = 0;
+   TString tempconst = "";
+   TString matname;
+   Bool_t composite = kFALSE;
+
+   if (z == 1) {
+      Double_t a = 0;
+      Double_t d = 0;
+
+      name = gdml->GetAttr(node, "name");
+      if ((strcmp(fCurrentFile, fStartFile)) != 0) {
+         name = TString::Format("%s_%s", name.Data(), fCurrentFile);
+      }
 
-  TGeoManager* mgr = gGeoManager;
-  TGeoElementTable* tab_ele = mgr->GetElementTable();
-  // We have to assume the media are monotonic increasing starting with 1
-  static int medid = mgr->GetListOfMedia()->GetSize()+1;
-  XMLNodePointer_t child = gdml->GetChild(node);
-  TString tempattr = "";
-  Int_t ncompo = 0, mixflag = 2;
-  Double_t density = 0;
-  TString name = "";
-  TGeoMixture* mix = 0;
-  TGeoMaterial* mat = 0;
-  TString tempconst = "";
-  TString matname;
-  Bool_t composite = kFALSE;
+      while (child != 0) {
+         attr = gdml->GetFirstAttr(child);
 
-  if (z == 1) {
-    Double_t a = 0;
-    Double_t d = 0;
+         if ((strcmp(gdml->GetNodeName(child), "property")) == 0) {
+            TNamed *property = new TNamed();
+            while (attr != 0) {
+               tempattr = gdml->GetAttrName(attr);
+               tempattr.ToLower();
 
-    while (child != 0) {
-      attr = gdml->GetFirstAttr(child);
+               if (tempattr == "name") {
+                  property->SetName(gdml->GetAttrValue(attr));
+               }
+               else if(tempattr == "ref") {
+                  property->SetTitle(gdml->GetAttrValue(attr));
+                  TGDMLMatrix *matrix = fmatrices[property->GetTitle()];
+                  if (!matrix)
+                     Error("MatProcess", "Reference matrix %s for material %s not found", property->GetTitle(), name.Data());
+                  properties.Add(property);
+               }
+               attr = gdml->GetNextAttr(attr);
+            }
+         }
 
-      if ((strcmp(gdml->GetNodeName(child), "atom")) == 0) {
-        while (attr != 0) {
-          tempattr = gdml->GetAttrName(attr);
-          tempattr.ToLower();
+         if ((strcmp(gdml->GetNodeName(child), "atom")) == 0) {
+            while (attr != 0) {
+               tempattr = gdml->GetAttrName(attr);
+               tempattr.ToLower();
 
-          if (tempattr == "value") {
-            a = Value(gdml->GetAttrValue(attr));
-          }
-          attr = gdml->GetNextAttr(attr);
-        }
-      }
+               if (tempattr == "value") {
+                  a = Value(gdml->GetAttrValue(attr));
+               }
+               attr = gdml->GetNextAttr(attr);
+            }
+         }
 
-      if ((strcmp(gdml->GetNodeName(child), "D")) == 0) {
-        while (attr != 0) {
-          tempattr = gdml->GetAttrName(attr);
-          tempattr.ToLower();
+         if ((strcmp(gdml->GetNodeName(child), "D")) == 0) {
+            while (attr != 0) {
+               tempattr = gdml->GetAttrName(attr);
+               tempattr.ToLower();
 
-          if (tempattr == "value") {
-            d = Value(gdml->GetAttrValue(attr));
-          }
-          attr = gdml->GetNextAttr(attr);
-        }
+               if (tempattr == "value") {
+                  d = Value(gdml->GetAttrValue(attr));
+               }
+               attr = gdml->GetNextAttr(attr);
+            }
+         }
+         child = gdml->GetNext(child);
       }
-      child = gdml->GetNext(child);
-    }
-    //still in the is Z else...but not in the while..
+      //still in the is Z else...but not in the while..
+      //CHECK FOR CONSTANTS
+      tempconst = gdml->GetAttr(node, "Z");
 
-    name = gdml->GetAttr(node, "name");
+      Double_t valZ = Value(tempconst);
 
-    if ((strcmp(fCurrentFile, fStartFile)) != 0) {
-      name = TString::Format("%s_%s", name.Data(), fCurrentFile);
-    }
-
-    //CHECK FOR CONSTANTS
-    tempconst = gdml->GetAttr(node, "Z");
+      TString tmpname = name;
+      //deal with special case - Z of vacuum is always 0
+      tmpname.ToLower();
+      if (tmpname == "vacuum") {
+         valZ = 0;
+      }
+      TString mat_name = NameShort(name);
+      mat = mgr->GetMaterial(mat_name);
+      if ( !mat )  {
+         mat = new TGeoMaterial(mat_name, a, valZ, d);
+      }
+      else  {
+         Info("TGDMLParse","Re-use existing material: %s",mat->GetName());
+      }
+      if (properties.GetSize()) {
+         TNamed *property;
+         TIter next(&properties);
+         while ((property = (TNamed*)next()))
+            mat->AddProperty(property->GetName(), property->GetTitle());
+      }
+      mixflag = 0;
+      //Note: Object(name, title) corresponds to Element(formula, name)
+      TGeoElement* mat_ele = tab_ele->FindElement(mat_name);
+      // We cannot use elements with Z = 0, so we expect a user definition
+      if (mat_ele && mat_ele->Z() == 0)
+         mat_ele = nullptr;
 
-    Double_t valZ = Value(tempconst);
+      if ( !mat_ele )  {
+         mat_ele = new TGeoElement(mat_name, mat_name, atoi(tempconst), a);
+      }
+      else if ( gDebug >= 2 )  {
+         Info("TGDMLParse","Re-use existing material-element: %s",mat_ele->GetName());
+      }
+      felemap[name.Data()] = mat_ele;
+   }
 
-    TString tmpname = name;
-    //deal with special case - Z of vacuum is always 0
-    tmpname.ToLower();
-    if (tmpname == "vacuum") {
-      valZ = 0;
-    }
-    TString mat_name = NameShort(name);
-    mat = mgr->GetMaterial(mat_name);
-    if ( !mat )  {
-      mat = new TGeoMaterial(mat_name, a, valZ, d);
-    }
-    else  {
-      Info("TGDMLParse","Re-use existing material: %s",mat->GetName());
-    }
-    mixflag = 0;
-    //Note: Object(name, title) corresponds to Element(formula, name)
-    TGeoElement* mat_ele = tab_ele->FindElement(mat_name);
-    // We cannot use elements with Z = 0, so we expect a user definition
-    if (mat_ele && mat_ele->Z() == 0)
-       mat_ele = nullptr;
+   else if (z == 0) {
+      while (child != 0) {
+         attr = gdml->GetFirstAttr(child);
 
-    if ( !mat_ele )  {
-      mat_ele = new TGeoElement(mat_name, mat_name, atoi(tempconst), a);
-    }
-    else if ( gDebug >= 2 )  {
-       Info("TGDMLParse","Re-use existing material-element: %s",mat_ele->GetName());
-    }
-    felemap[name.Data()] = mat_ele;
-  }
+         if ((strcmp(gdml->GetNodeName(child), "property")) == 0) {
+            TNamed *property = new TNamed();
+            while (attr != 0) {
+               tempattr = gdml->GetAttrName(attr);
+               tempattr.ToLower();
 
-  else if (z == 0) {
-     while (child != 0) {
-        attr = gdml->GetFirstAttr(child);
+               if (tempattr == "name") {
+                  property->SetName(gdml->GetAttrValue(attr));
+               }
+               else if(tempattr == "ref") {
+                  property->SetTitle(gdml->GetAttrValue(attr));
+                  TGDMLMatrix *matrix = fmatrices[property->GetTitle()];
+                  if (!matrix)
+                     Error("MatProcess", "Reference matrix %s for material %s not found", property->GetTitle(), name.Data());
+                  properties.Add(property);
+               }
+               attr = gdml->GetNextAttr(attr);
+            }
+         }
+         if ((strcmp(gdml->GetNodeName(child), "fraction")) == 0) {
+            Double_t n = 0;
+            TString ref = "";
+            ncompo = ncompo + 1;
 
-        if ((strcmp(gdml->GetNodeName(child), "fraction")) == 0) {
-           Double_t n = 0;
-           TString ref = "";
-           ncompo = ncompo + 1;
+            while (attr != 0) {
+               tempattr = gdml->GetAttrName(attr);
+               tempattr.ToLower();
 
-           while (attr != 0) {
-              tempattr = gdml->GetAttrName(attr);
-              tempattr.ToLower();
+               if (tempattr == "n") {
+                  n = Value(gdml->GetAttrValue(attr));
+               } else if (tempattr == "ref") {
+                  ref = gdml->GetAttrValue(attr);
+                  if ((strcmp(fCurrentFile, fStartFile)) != 0) {
+                     ref = TString::Format("%s_%s", ref.Data(), fCurrentFile);
+                  }
+               }
+               attr = gdml->GetNextAttr(attr);
+            }
+            fracmap[ref.Data()] = n;
+         }
 
-              if (tempattr == "n") {
-                 n = Value(gdml->GetAttrValue(attr));
-              } else if (tempattr == "ref") {
-                 ref = gdml->GetAttrValue(attr);
-                 if ((strcmp(fCurrentFile, fStartFile)) != 0) {
-                    ref = TString::Format("%s_%s", ref.Data(), fCurrentFile);
-                 }
+         else if ((strcmp(gdml->GetNodeName(child), "composite")) == 0) {
+            composite = kTRUE;
+            Double_t n = 0;
+            TString ref = "";
+            ncompo = ncompo + 1;
 
-              }
+            while (attr != 0) {
+               tempattr = gdml->GetAttrName(attr);
+               tempattr.ToLower();
+               if (tempattr == "n") {
+                  n = Value(gdml->GetAttrValue(attr));
+               } else if (tempattr == "ref") {
+                  ref = gdml->GetAttrValue(attr);
+                  if ((strcmp(fCurrentFile, fStartFile)) != 0) {
+                     ref = TString::Format("%s_%s", ref.Data(), fCurrentFile);
+                  }
+               }
+               attr = gdml->GetNextAttr(attr);
+            }
+            fracmap[ref.Data()] = n;
+         }
+         else if ((strcmp(gdml->GetNodeName(child), "D")) == 0) {
+            while (attr != 0) {
+               tempattr = gdml->GetAttrName(attr);
+               tempattr.ToLower();
 
-              attr = gdml->GetNextAttr(attr);
-           }
-           fracmap[ref.Data()] = n;
+               if (tempattr == "value") {
+                  density = Value(gdml->GetAttrValue(attr));
+               }
+               attr = gdml->GetNextAttr(attr);
+            }
+         }
+         child = gdml->GetNext(child);
+      }
+      //still in the not Z else...but not in the while..
 
-        }
+      name = gdml->GetAttr(node, "name");
+      if ((strcmp(fCurrentFile, fStartFile)) != 0) {
+         name = TString::Format("%s_%s", name.Data(), fCurrentFile);
+      }
+      //mix = new TGeoMixture(NameShort(name), 0 /*ncompo*/, density);
+      mixflag = 1;
+      TString mat_name = NameShort(name);
+      mat = mgr->GetMaterial(mat_name);
+      if ( !mat )  {
+         mix = new TGeoMixture(mat_name, ncompo, density);
+      }
+      else if ( mat->IsMixture() ) {
+         mix = (TGeoMixture*)mat;
+         if ( gDebug >= 2 )
+            Info("TGDMLParse","Re-use existing material-mixture: %s",mix->GetName());
+      }
+      else  {
+         Error("TGDMLParse","WARNING! Inconsistent material definitions between GDML and TGeoManager");
+      }
+      if (properties.GetSize()) {
+         TNamed *property;
+         TIter next(&properties);
+         while ((property = (TNamed*)next()))
+            mix->AddProperty(property->GetName(), property->GetTitle());
+      }
+      Int_t natoms;
+      Double_t weight;
 
-        else if ((strcmp(gdml->GetNodeName(child), "composite")) == 0) {
-           composite = kTRUE;
-           Double_t n = 0;
-           TString ref = "";
-           ncompo = ncompo + 1;
-
-           while (attr != 0) {
-              tempattr = gdml->GetAttrName(attr);
-              tempattr.ToLower();
-
-              if (tempattr == "n") {
-                 n = Value(gdml->GetAttrValue(attr));
-              } else if (tempattr == "ref") {
-                 ref = gdml->GetAttrValue(attr);
-                 if ((strcmp(fCurrentFile, fStartFile)) != 0) {
-                    ref = TString::Format("%s_%s", ref.Data(), fCurrentFile);
-                 }
-              }
-
-              attr = gdml->GetNextAttr(attr);
-           }
+      for (fractions f = fracmap.begin(); f != fracmap.end(); ++f) {
+         matname = f->first;
+         matname = NameShort(matname);
 
-           fracmap[ref.Data()] = n;
+         TGeoMaterial *mattmp = (TGeoMaterial*)gGeoManager->GetListOfMaterials()->FindObject(matname);
 
-        }
-        else if ((strcmp(gdml->GetNodeName(child), "D")) == 0) {
-           while (attr != 0) {
-              tempattr = gdml->GetAttrName(attr);
-              tempattr.ToLower();
+         if (mattmp || (felemap.find(f->first) != felemap.end())) {
+            if (composite) {
+               natoms = (Int_t)f->second;
 
-              if (tempattr == "value") {
-                 density = Value(gdml->GetAttrValue(attr));
-              }
+               mix->AddElement(felemap[f->first], natoms);
 
-              attr = gdml->GetNextAttr(attr);
-           }
-        }
+            }
 
-        child = gdml->GetNext(child);
-     }
-     //still in the not Z else...but not in the while..
+            else {
+               weight = f->second;
+               if (mattmp){
+                  mix->AddElement(mattmp, weight);
+               }
+               else {
+                  mix->AddElement(felemap[f->first], weight);
+               }
+            }
+         }
+      }
+   }//end of not Z else
 
-     name = gdml->GetAttr(node, "name");
-     if ((strcmp(fCurrentFile, fStartFile)) != 0) {
-        name = TString::Format("%s_%s", name.Data(), fCurrentFile);
-     }
-     //mix = new TGeoMixture(NameShort(name), 0 /*ncompo*/, density);
-     mixflag = 1;
-     TString mat_name = NameShort(name);
-     mat = mgr->GetMaterial(mat_name);
-     if ( !mat )  {
-       mix = new TGeoMixture(mat_name, ncompo, density);
-     }
-     else if ( mat->IsMixture() ) {
-       mix = (TGeoMixture*)mat;
-       if ( gDebug >= 2 )
-          Info("TGDMLParse","Re-use existing material-mixture: %s",mix->GetName());
-     }
-     else  {
-       Error("TGDMLParse","WARNING! Inconsistent material definitions between GDML and TGeoManager");
-     }
-     Int_t natoms;
-     Double_t weight;
+   medid = medid + 1;
 
-     for (fractions f = fracmap.begin(); f != fracmap.end(); ++f) {
-        matname = f->first;
-        matname = NameShort(matname);
+   TGeoMedium* med = mgr->GetMedium(NameShort(name));
+   if ( !med )   {
+      if (mixflag == 1) {
+         fmixmap[name.Data()] = mix;
+         med = new TGeoMedium(NameShort(name), medid, mix);
+      } else if (mixflag == 0) {
+         fmatmap[name.Data()] = mat;
+         med = new TGeoMedium(NameShort(name), medid, mat);
+      }
+   }
+   else if ( gDebug >= 2 ) {
+      Info("TGDMLParse","Re-use existing medium: %s",med->GetName());
+   }
+   fmedmap[name.Data()] = med;
 
-        TGeoMaterial *mattmp = (TGeoMaterial*)gGeoManager->GetListOfMaterials()->FindObject(matname);
+   return child;
+}
 
-        if (mattmp || (felemap.find(f->first) != felemap.end())) {
-           if (composite) {
-              natoms = (Int_t)f->second;
+////////////////////////////////////////////////////////////////////////////////
+/// In the structure section of the GDML file, skin surfaces can be declared.
 
-              mix->AddElement(felemap[f->first], natoms);
+XMLNodePointer_t TGDMLParse::SkinSurfaceProcess(TXMLEngine* gdml, XMLNodePointer_t node, XMLAttrPointer_t attr)
+{
+   TString name, surfname, volname;
+   TString tempattr;
 
-           }
+   while (attr != 0) {
+      tempattr = gdml->GetAttrName(attr);
+      tempattr.ToLower();
 
-           else {
-              weight = f->second;
-              if (mattmp){
+      if (tempattr == "name") {
+         name = gdml->GetAttrValue(attr);
+      }
+      if (tempattr == "surfaceproperty") {
+         surfname = gdml->GetAttrValue(attr);
+      }
+      attr = gdml->GetNextAttr(attr);
+   }
 
-                 mix->AddElement(mattmp, weight);
-              }
-              else {
+   XMLNodePointer_t child = gdml->GetChild(node);
+   while (child != 0) {
+      attr = gdml->GetFirstAttr(child);
+      if ((strcmp(gdml->GetNodeName(child), "volumeref")) == 0) {
+         while (attr != 0) {
+            tempattr = gdml->GetAttrName(attr);
+            tempattr.ToLower();
+            if (tempattr == "ref") {
+               volname = gdml->GetAttrValue(attr);
+            }
+            attr = gdml->GetNextAttr(attr);
+         }
+      } // loop on child attributes
+      child = gdml->GetNext(child);
+   } // loop on children
+   TGeoOpticalSurface *surf = gGeoManager->GetOpticalSurface(surfname);
+   if (!surf) 
+      Fatal("SkinSurfaceProcess", "Skin surface %s: referenced optical surface %s not defined",
+            name.Data(), surfname.Data());
+   TGeoVolume *vol = fvolmap[volname.Data()];
+   TGeoSkinSurface *skin = new TGeoSkinSurface(name, surfname, surf, vol);
+   gGeoManager->AddSkinSurface(skin);
+   return child;
+}
 
-                 mix->AddElement(felemap[f->first], weight);
-              }
-           }
-        }
-     }
+////////////////////////////////////////////////////////////////////////////////
+/// In the structure section of the GDML file, border surfaces can be declared.
 
-  }//end of not Z else
+XMLNodePointer_t TGDMLParse::BorderSurfaceProcess(TXMLEngine* gdml, XMLNodePointer_t node, XMLAttrPointer_t attr)
+{
+   TString name, surfname, nodename[2];
+   TString tempattr;
 
-   medid = medid + 1;
+   while (attr != 0) {
+      tempattr = gdml->GetAttrName(attr);
+      tempattr.ToLower();
 
-   TGeoMedium* med = mgr->GetMedium(NameShort(name));
-   if ( !med )   {
-     if (mixflag == 1) {
-       fmixmap[name.Data()] = mix;
-       med = new TGeoMedium(NameShort(name), medid, mix);
-     } else if (mixflag == 0) {
-       fmatmap[name.Data()] = mat;
-       med = new TGeoMedium(NameShort(name), medid, mat);
-     }
-   }
-   else if ( gDebug >= 2 ) {
-      Info("TGDMLParse","Re-use existing medium: %s",med->GetName());
+      if (tempattr == "name") {
+         name = gdml->GetAttrValue(attr);
+      }
+      if (tempattr == "surfaceproperty") {
+         surfname = gdml->GetAttrValue(attr);
+      }
+      attr = gdml->GetNextAttr(attr);
    }
-   fmedmap[name.Data()] = med;
 
+   XMLNodePointer_t child = gdml->GetChild(node);
+   Int_t inode = 0;
+   while (child != 0) {
+      attr = gdml->GetFirstAttr(child);
+      if ((strcmp(gdml->GetNodeName(child), "physvolref")) == 0) {
+         while (attr != 0) {
+            tempattr = gdml->GetAttrName(attr);
+            tempattr.ToLower();
+            if (tempattr == "ref") {
+               nodename[inode++] = gdml->GetAttrValue(attr);
+            }
+            attr = gdml->GetNextAttr(attr);
+         }
+      } // loop on child attributes
+      child = gdml->GetNext(child);
+   } // loop on children
+   if (inode != 2)
+      Fatal("BorderSurfaceProcess", "Border surface %s not referencing two nodes", name.Data());
+   TGeoOpticalSurface *surf = gGeoManager->GetOpticalSurface(surfname);
+   if (!surf)
+      Fatal("BorderSurfaceProcess", "Border surface %s: referenced optical surface %s not defined",
+            name.Data(), surfname.Data());
+   TGeoNode *node1 = fpvolmap[nodename[0].Data()];
+   TGeoNode *node2 = fpvolmap[nodename[1].Data()];
+   if (!node1 || !node2)
+      Fatal("BorderSurfaceProcess", "Border surface %s: not found nodes %s or %s",
+            name.Data(), nodename[0].Data(), nodename[1].Data());
+
+   TGeoBorderSurface *border = new TGeoBorderSurface(name, surfname, surf, node1, node2);
+   gGeoManager->AddBorderSurface(border);
    return child;
-
 }
 
 ////////////////////////////////////////////////////////////////////////////////
@@ -1562,8 +1866,10 @@ XMLNodePointer_t TGDMLParse::VolProcess(TXMLEngine* gdml, XMLNodePointer_t node)
 // END: reflectedSolid
 
          vol->AddNode(lv, copynum, transform);
+         TGeoNode *lastnode = (TGeoNode*)vol->GetNodes()->Last();
          if (!pnodename.IsNull())
-            ((TNamed*)vol->GetNodes()->Last())->SetName(pnodename);
+            lastnode->SetName(pnodename);
+         fpvolmap[lastnode->GetName()] = lastnode;
       } else if ((strcmp(gdml->GetNodeName(child), "divisionvol")) == 0) {
 
          TString divVolref = "";
@@ -2108,9 +2414,10 @@ XMLNodePointer_t TGDMLParse::AssProcess(TXMLEngine* gdml, XMLNodePointer_t node)
          fVolID = fVolID + 1;
          matr = new TGeoCombiTrans(*pos, *rot);
          assem->AddNode(lv, copynum, matr);
+         TGeoNode *lastnode = (TGeoNode*)assem->GetNodes()->Last();
          if (!pnodename.IsNull())
-            ((TNamed*)assem->GetNodes()->Last())->SetName(pnodename);
-
+            lastnode->SetName(pnodename);
+         fpvolmap[lastnode->GetName()] = lastnode;
       }
       child = gdml->GetNext(child);
    }
diff --git a/geom/gdml/src/TGDMLWrite.cxx b/geom/gdml/src/TGDMLWrite.cxx
index c64a1af3a75..de095dcb62c 100644
--- a/geom/gdml/src/TGDMLWrite.cxx
+++ b/geom/gdml/src/TGDMLWrite.cxx
@@ -126,6 +126,8 @@ See that function for details.
 
 */
 
+#include "TGDMLWrite.h"
+
 #include "TGeoManager.h"
 #include "TGeoMaterial.h"
 #include "TGeoMatrix.h"
@@ -155,7 +157,7 @@ See that function for details.
 #include "TGeoElement.h"
 #include "TGeoShape.h"
 #include "TGeoCompositeShape.h"
-#include "TGDMLWrite.h"
+#include "TGeoOpticalSurface.h"
 #include <stdlib.h>
 #include <string>
 #include <map>
@@ -351,6 +353,7 @@ void TGDMLWrite::WriteGDMLfile(TGeoManager * geomanager,
    //calling main extraction functions (with measuring time)
    time_t startT, endT;
    startT = time(NULL);
+   ExtractMatrices(geomanager->GetListOfGDMLMatrices());
    fMaterialsNode = ExtractMaterials(materialsLst);
 
    Info("WriteGDMLfile", "Extracting volumes");
@@ -358,6 +361,9 @@ void TGDMLWrite::WriteGDMLfile(TGeoManager * geomanager,
    Info("WriteGDMLfile", "%i solids added", fSolCnt);
    Info("WriteGDMLfile", "%i volumes added", fVolCnt);
    Info("WriteGDMLfile", "%i physvolumes added", fPhysVolCnt);
+   ExtractOpticalSurfaces(geomanager->GetListOfOpticalSurfaces());
+   ExtractSkinSurfaces(geomanager->GetListOfSkinSurfaces());
+   ExtractBorderSurfaces(geomanager->GetListOfBorderSurfaces());
    endT = time(NULL);
    //<gdml>
    fGdmlE->AddChild(rootNode, fDefineNode);                 //  <define>...</define>
@@ -381,6 +387,65 @@ void TGDMLWrite::WriteGDMLfile(TGeoManager * geomanager,
    delete fGdmlE;
 }
 
+////////////////////////////////////////////////////////////////////////////////
+/// Method exporting GDML matrices
+
+void TGDMLWrite::ExtractMatrices(TObjArray* matrixList)
+{
+   if (!matrixList->GetEntriesFast()) return;
+   XMLNodePointer_t matrixN;
+   TIter next(matrixList);
+   TGDMLMatrix *matrix;
+   while ((matrix = (TGDMLMatrix*)next())) {
+      matrixN = CreateMatrixN(matrix);
+      fGdmlE->AddChild(fDefineNode, matrixN);
+   }
+}
+
+////////////////////////////////////////////////////////////////////////////////
+/// Method exporting optical surfaces
+
+void TGDMLWrite::ExtractOpticalSurfaces(TObjArray *surfaces)
+{
+   if (!surfaces->GetEntriesFast()) return;
+   XMLNodePointer_t surfaceN;
+   TIter next(surfaces);
+   TGeoOpticalSurface *surf;
+   while ((surf = (TGeoOpticalSurface*)next())) {
+      surfaceN = CreateOpticalSurfaceN(surf);
+      fGdmlE->AddChild(fSolidsNode, surfaceN);
+   }
+}
+
+////////////////////////////////////////////////////////////////////////////////
+/// Method exporting skin surfaces
+
+void TGDMLWrite::ExtractSkinSurfaces(TObjArray *surfaces)
+{
+   if (!surfaces->GetEntriesFast()) return;
+   XMLNodePointer_t surfaceN;
+   TIter next(surfaces);
+   TGeoSkinSurface *surf;
+   while ((surf = (TGeoSkinSurface*)next())) {
+      surfaceN = CreateSkinSurfaceN(surf);
+      fGdmlE->AddChild(fStructureNode, surfaceN);
+   }
+}
+
+////////////////////////////////////////////////////////////////////////////////
+/// Method exporting border surfaces
+
+void TGDMLWrite::ExtractBorderSurfaces(TObjArray *surfaces)
+{
+   if (!surfaces->GetEntriesFast()) return;
+   XMLNodePointer_t surfaceN;
+   TIter next(surfaces);
+   TGeoBorderSurface *surf;
+   while ((surf = (TGeoBorderSurface*)next())) {
+      surfaceN = CreateBorderSurfaceN(surf);
+      fGdmlE->AddChild(fStructureNode, surfaceN);
+   }
+}
 
 ////////////////////////////////////////////////////////////////////////////////
 /// Method exporting materials
@@ -629,6 +694,17 @@ XMLNodePointer_t TGDMLWrite::CreateFractionN(Double_t percentage, const char * r
    return fractN;
 }
 
+////////////////////////////////////////////////////////////////////////////////
+/// Creates "property" node for GDML
+
+XMLNodePointer_t TGDMLWrite::CreatePropertyN(TNamed const &property)
+{
+  XMLNodePointer_t propertyN = fGdmlE->NewChild(0, 0, "property", 0);
+  fGdmlE->NewAttr(propertyN, 0, "name", property.GetName());
+  fGdmlE->NewAttr(propertyN, 0, "ref", property.GetTitle());
+  return propertyN;
+}
+
 ////////////////////////////////////////////////////////////////////////////////
 /// Creates "isotope" node for GDML
 
@@ -756,6 +832,15 @@ XMLNodePointer_t TGDMLWrite::CreateMixtureN(TGeoMixture * mixture, XMLNodePointe
       fGdmlE->AddChild(mainN, CreateFractionN(wPercentage[itr->first], itr->first.Data()));
    }
 
+   // Write properties
+   TList const &properties = mixture->GetProperties();
+   if (properties.GetSize()) {
+      TIter next(&properties);
+      TNamed *property;
+      while ((property = (TNamed*)next()))
+        fGdmlE->AddChild(mainN, CreatePropertyN(*property));
+   }
+
    return mainN;
 }
 
@@ -787,6 +872,14 @@ XMLNodePointer_t TGDMLWrite::CreateMaterialN(TGeoMaterial * material, TString mn
    fGdmlE->NewAttr(mainN, 0, "Z", TString::Format(fltPrecision.Data(), valZ)); //material->GetZ()));
    fGdmlE->AddChild(mainN, CreateDN(material->GetDensity()));
    fGdmlE->AddChild(mainN, CreateAtomN(material->GetA()));
+   // Create properties if any
+   TList const &properties = material->GetProperties();
+   if (properties.GetSize()) {
+      TIter next(&properties);
+      TNamed *property;
+      while ((property = (TNamed*)next()))
+        fGdmlE->AddChild(mainN, CreatePropertyN(*property));
+   }
    return mainN;
 }
 
@@ -1613,6 +1706,61 @@ XMLNodePointer_t TGDMLWrite::CreateCommonBoolN(TGeoCompositeShape *geoShape)
    return mainN;
 }
 
+////////////////////////////////////////////////////////////////////////////////
+/// Creates "opticalsurface" node for GDML
+
+XMLNodePointer_t TGDMLWrite::CreateOpticalSurfaceN(TGeoOpticalSurface * geoSurf)
+{
+   XMLNodePointer_t mainN = fGdmlE->NewChild(0, 0, "opticalsurface", 0);
+   const TString fltPrecision = TString::Format("%%.%dg", fFltPrecision);
+   fGdmlE->NewAttr(mainN, 0, "name", geoSurf->GetName());
+   fGdmlE->NewAttr(mainN, 0, "model", TGeoOpticalSurface::ModelToString(geoSurf->GetModel()));
+   fGdmlE->NewAttr(mainN, 0, "finish", TGeoOpticalSurface::FinishToString(geoSurf->GetFinish()));
+   fGdmlE->NewAttr(mainN, 0, "type", TGeoOpticalSurface::TypeToString(geoSurf->GetType()));
+   fGdmlE->NewAttr(mainN, 0, "value", TString::Format(fltPrecision.Data(), geoSurf->GetValue()));
+
+   // Write properties
+   TList const &properties = geoSurf->GetProperties();
+   if (properties.GetSize()) {
+      TIter next(&properties);
+      TNamed *property;
+      while ((property = (TNamed*)next()))
+        fGdmlE->AddChild(mainN, CreatePropertyN(*property));
+   }
+
+   return mainN;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+/// Creates "skinsurface" node for GDML
+
+XMLNodePointer_t TGDMLWrite::CreateSkinSurfaceN(TGeoSkinSurface * geoSurf)
+{
+   XMLNodePointer_t mainN = fGdmlE->NewChild(0, 0, "skinsurface", 0);
+   fGdmlE->NewAttr(mainN, 0, "name", geoSurf->GetName());
+   fGdmlE->NewAttr(mainN, 0, "surfaceproperty", geoSurf->GetTitle());
+   // Cretate the logical volume reference node
+   auto childN = fGdmlE->NewChild(0, 0, "volumeref", 0);
+   fGdmlE->NewAttr(childN, 0, "ref", geoSurf->GetVolume()->GetName());
+   fGdmlE->AddChild(mainN, childN);
+   return mainN;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+/// Creates "bordersurface" node for GDML
+
+XMLNodePointer_t TGDMLWrite::CreateBorderSurfaceN(TGeoBorderSurface * geoSurf)
+{
+   XMLNodePointer_t mainN = fGdmlE->NewChild(0, 0, "bordersurface", 0);
+   fGdmlE->NewAttr(mainN, 0, "name", geoSurf->GetName());
+   fGdmlE->NewAttr(mainN, 0, "surfaceproperty", geoSurf->GetTitle());
+   // Cretate the logical volume reference node
+   auto childN = fGdmlE->NewChild(0, 0, "physvolref", 0);
+   fGdmlE->NewAttr(childN, 0, "ref", geoSurf->GetNode1()->GetName());
+   fGdmlE->NewAttr(childN, 0, "ref", geoSurf->GetNode2()->GetName());
+   fGdmlE->AddChild(mainN, childN);
+   return mainN;
+}
 
 ////////////////////////////////////////////////////////////////////////////////
 /// Creates "position" kind of node for GDML
@@ -1644,6 +1792,18 @@ XMLNodePointer_t TGDMLWrite::CreateRotationN(const char * name, Xyz rotation, co
    return mainN;
 }
 
+////////////////////////////////////////////////////////////////////////////////
+/// Creates "matrix" kind of node for GDML
+
+XMLNodePointer_t TGDMLWrite::CreateMatrixN(TGDMLMatrix const *matrix)
+{
+   XMLNodePointer_t mainN = fGdmlE->NewChild(0, 0, "matrix", 0);
+   fGdmlE->NewAttr(mainN, 0, "name", matrix->GetName());
+   fGdmlE->NewAttr(mainN, 0, "coldim", TString::Format("%ld", matrix->GetCols()));
+   fGdmlE->NewAttr(mainN, 0, "values", matrix->GetMatrixAsString());
+   return mainN;
+}
+
 ////////////////////////////////////////////////////////////////////////////////
 /// Creates "setup" node for GDML
 
diff --git a/geom/geom/CMakeLists.txt b/geom/geom/CMakeLists.txt
index 520be8b5d59..6bcc26028ae 100644
--- a/geom/geom/CMakeLists.txt
+++ b/geom/geom/CMakeLists.txt
@@ -3,7 +3,7 @@
 ############################################################################
 
 set(headers1 TGeoAtt.h TGeoStateInfo.h TGeoBoolNode.h
-             TGeoMedium.h TGeoMaterial.h
+             TGeoMedium.h TGeoMaterial.h TGeoOpticalSurface.h
              TGeoMatrix.h TGeoVolume.h TGeoNode.h
              TGeoVoxelFinder.h TGeoShape.h TGeoBBox.h
              TGeoPara.h TGeoTube.h TGeoTorus.h TGeoSphere.h
@@ -14,7 +14,7 @@ set(headers1 TGeoAtt.h TGeoStateInfo.h TGeoBoolNode.h
              TVirtualGeoConverter.h TGeoPolygon.h TGeoXtru.h TGeoPhysicalNode.h
              TGeoHelix.h TGeoParaboloid.h TGeoElement.h TGeoHalfSpace.h
              TGeoBuilder.h TGeoNavigator.h TGeoRegion.h
-             TGeoPhysicalConstants.h TGeoSystemOfUnits.h)
+             TGeoPhysicalConstants.h TGeoSystemOfUnits.h TGDMLMatrix.h)
 set(headers2 TGeoPatternFinder.h TGeoCache.h TVirtualMagField.h
              TGeoUniformMagField.h TGeoGlobalMagField.h TGeoBranchArray.h
              TGeoExtension.h TGeoParallelWorld.h)
diff --git a/geom/geom/inc/LinkDef1.h b/geom/geom/inc/LinkDef1.h
index f8b237d7140..eb0655fde41 100644
--- a/geom/geom/inc/LinkDef1.h
+++ b/geom/geom/inc/LinkDef1.h
@@ -20,6 +20,12 @@
 #pragma link C++ class TGeoIntersection+;
 #pragma link C++ class TGeoSubtraction+;
 #pragma link C++ class TGeoMedium+;
+#pragma link C++ class TGeoOpticalSurface+;
+#pragma link C++ enum  TGeoOpticalSurface::ESurfaceType;
+#pragma link C++ enum  TGeoOpticalSurface::ESurfaceModel;
+#pragma link C++ enum  TGeoOpticalSurface::ESurfaceFinish;
+#pragma link C++ class TGeoSkinSurface+;
+#pragma link C++ class TGeoBorderSurface+;
 #pragma link C++ class TGeoElement+;
 #pragma read sourceClass="TGeoElement" targetClass="TGeoElement" version="[1-2]" source="" target="" \
     code="{ newObj->ComputeDerivedQuantities() ; }" 
@@ -90,6 +96,7 @@
 #pragma link C++ class TGeoBuilder;
 #pragma link C++ class TGeoNavigator+;
 #pragma link C++ class TGeoNavigatorArray;
+#pragma link C++ class TGDMLMatrix+;
 #pragma link C++ struct std::map<std::thread::id, TGeoNavigatorArray *>;
 #pragma link C++ struct std::pair<std::thread::id, TGeoNavigatorArray *>;
 #pragma link C++ struct std::map<std::thread::id, Int_t>;
diff --git a/geom/geom/inc/TGDMLMatrix.h b/geom/geom/inc/TGDMLMatrix.h
new file mode 100644
index 00000000000..a79d03db6b6
--- /dev/null
+++ b/geom/geom/inc/TGDMLMatrix.h
@@ -0,0 +1,61 @@
+// @(#)root/gdml:$Id$
+// Author: Andrei Gheata 05/12/2018
+
+/*************************************************************************
+ * Copyright (C) 1995-2011, Rene Brun and Fons Rademakers.               *
+ * All rights reserved.                                                  *
+ *                                                                       *
+ * For the licensing terms see $ROOTSYS/LICENSE.                         *
+ * For the list of contributors see $ROOTSYS/README/CREDITS.             *
+ *************************************************************************/
+
+#ifndef ROOT_TGDMLMATRIX
+#define ROOT_TGDMLMATRIX
+
+#include <TNamed.h>
+
+
+////////////////////////////////////////////////////////////////////////////
+//                                                                        //
+// TGDMLProperty - A property with a name and a reference name pointing   //
+//     to a GDML matrix object                                            //
+////////////////////////////////////////////////////////////////////////////
+
+typedef TNamed TGDMLProperty;
+
+////////////////////////////////////////////////////////////////////////////
+//                                                                        //
+// TGDMLMatrix - A matrix used for GDML parsing, the objects have to be   //
+//     exposed via TGeoManager interfcace to be able to construct optical //
+//     surfaces.                                                          //
+//                                                                        //
+////////////////////////////////////////////////////////////////////////////
+
+class TGDMLMatrix : public TNamed {
+public:
+   TGDMLMatrix() {}
+   TGDMLMatrix(const char *name, size_t rows,size_t cols);
+   TGDMLMatrix(const TGDMLMatrix& rhs);
+   TGDMLMatrix& operator=(const TGDMLMatrix& rhs);
+  ~TGDMLMatrix() { delete [] fMatrix; }
+
+   void        Set(size_t r, size_t c, Double_t a);
+   Double_t    Get(size_t r, size_t c) const;
+   size_t      GetRows() const { return fNrows; }
+   size_t      GetCols() const { return fNcols; }
+   void        SetMatrixAsString(const char *mat) { fTitle = mat; }
+   const char *GetMatrixAsString() const { return fTitle.Data(); }
+
+   void        Print(Option_t *option="") const;
+
+ private:
+
+   Int_t  fNelem = 0;                // Number of elements
+   size_t fNrows = 0;                // Number of rows
+   size_t fNcols = 0;                // Number of columns
+   Double_t *fMatrix = nullptr;      // [fNelem] Matrix elements
+
+   ClassDef(TGDMLMatrix, 1)          // Class representing a matrix used temporary for GDML parsing
+};
+
+#endif /* ROOT_TGDMLMATRIX */
diff --git a/geom/geom/inc/TGeoManager.h b/geom/geom/inc/TGeoManager.h
index 723081fd38f..25b9cf05021 100644
--- a/geom/geom/inc/TGeoManager.h
+++ b/geom/geom/inc/TGeoManager.h
@@ -34,6 +34,10 @@ class TVirtualGeoPainter;
 class THashList;
 class TGeoParallelWorld;
 class TGeoRegion;
+class TGDMLMatrix;
+class TGeoOpticalSurface;
+class TGeoSkinSurface;
+class TGeoBorderSurface;
 
 class TGeoManager : public TNamed
 {
@@ -89,7 +93,10 @@ private :
    TObjArray            *fGVolumes;         //! list of runtime volumes
    TObjArray            *fTracks;           //-> list of tracks attached to geometry
    TObjArray            *fPdgNames;         //-> list of pdg names for tracks
-//   TObjArray            *fNavigators;       //! list of navigators
+   TObjArray            *fGDMLMatrices;     //-> list of matrices read from GDML
+   TObjArray            *fOpticalSurfaces;  //-> list of optical surfaces read from GDML
+   TObjArray            *fSkinSurfaces;     //-> list of skin surfaces read from GDML
+   TObjArray            *fBorderSurfaces;   //-> list of border surfaces read from GDML
    TList                *fMaterials;        //-> list of materials
    TList                *fMedia;            //-> list of tracking media
    TObjArray            *fNodes;            //-> current branch of nodes
@@ -471,6 +478,10 @@ public:
    TObjArray             *GetListOfGShapes() const      {return fGShapes;}
    TObjArray             *GetListOfUVolumes() const     {return fUniqueVolumes;}
    TObjArray             *GetListOfTracks() const       {return fTracks;}
+   TObjArray             *GetListOfGDMLMatrices() const {return fGDMLMatrices;}
+   TObjArray             *GetListOfOpticalSurfaces() const {return fOpticalSurfaces;}
+   TObjArray             *GetListOfSkinSurfaces() const    {return fSkinSurfaces;}
+   TObjArray             *GetListOfBorderSurfaces() const  {return fBorderSurfaces;}
    TObjArray             *GetListOfRegions() const      {return fRegions;}
    TGeoNavigatorArray    *GetListOfNavigators() const;
    TGeoElementTable      *GetElementTable();
@@ -528,6 +539,17 @@ public:
    TGeoMedium            *GetMedium(const char *medium) const;
    TGeoMedium            *GetMedium(Int_t numed) const;
    Int_t                  GetMaterialIndex(const char *matname) const;
+
+   //--- GDML object accessors
+   TGDMLMatrix           *GetGDMLMatrix(const char *name) const;
+   void                   AddGDMLMatrix(TGDMLMatrix *mat);
+   TGeoOpticalSurface    *GetOpticalSurface(const char *name) const;
+   void                   AddOpticalSurface(TGeoOpticalSurface *optsurf);
+   TGeoSkinSurface       *GetSkinSurface(const char *name) const;
+   void                   AddSkinSurface(TGeoSkinSurface *surf);
+   TGeoBorderSurface     *GetBorderSurface(const char *name) const;
+   void                   AddBorderSurface(TGeoBorderSurface *surf);
+
 //   TGeoShape             *GetShape(const char *name) const;
    TGeoVolume            *GetVolume(const char*name) const;
    TGeoVolume            *GetVolume(Int_t uid) const {return (TGeoVolume*)fUniqueVolumes->At(uid);}
@@ -556,7 +578,7 @@ public:
    void                  SetUseParallelWorldNav(Bool_t flag);
    Bool_t                IsParallelWorldNav() const {return fUsePWNav;}
 
-   ClassDef(TGeoManager, 15)          // geometry manager
+   ClassDef(TGeoManager, 16)          // geometry manager
 };
 
 R__EXTERN TGeoManager *gGeoManager;
diff --git a/geom/geom/inc/TGeoMaterial.h b/geom/geom/inc/TGeoMaterial.h
index f554e533ca9..5325832f10c 100644
--- a/geom/geom/inc/TGeoMaterial.h
+++ b/geom/geom/inc/TGeoMaterial.h
@@ -12,13 +12,12 @@
 #ifndef ROOT_TGeoMaterial
 #define ROOT_TGeoMaterial
 
-#include "TNamed.h"
-
-#include "TAttFill.h"
+#include <TNamed.h>
+#include <TAttFill.h>
+#include <TList.h>
 
 #include "TGeoElement.h"
 
-
 // forward declarations
 class TGeoExtension;
 
@@ -54,6 +53,7 @@ protected:
    TObject                 *fShader;     // shader with optical properties
    TObject                 *fCerenkov;   // pointer to class with Cerenkov properties
    TGeoElement             *fElement;    // pointer to element composing the material
+   TList                    fProperties; // user-defined properties
    TGeoExtension           *fUserExtension;  //! Transient user-defined extension to materials
    TGeoExtension           *fFWExtension;    //! Transient framework-defined extension to materials
 
@@ -80,6 +80,9 @@ public:
    virtual TGeoMaterial    *DecayMaterial(Double_t time, Double_t precision=0.001);
    virtual void             FillMaterialEvolution(TObjArray *population, Double_t precision=0.001);
    // getters & setters
+   bool                     AddProperty(const char *property, const char *ref);
+   const char              *GetPropertyRef(const char *property);
+   TList const             &GetProperties() const { return fProperties; }
    virtual Int_t            GetByteCount() const {return sizeof(*this);}
    virtual Double_t         GetA() const       {return fA;}
    virtual Double_t         GetZ()  const      {return fZ;}
@@ -125,7 +128,7 @@ public:
 
 
 
-   ClassDef(TGeoMaterial, 5)              // base material class
+   ClassDef(TGeoMaterial, 6)              // base material class
 
 //***** Need to add classes and globals to LinkDef.h *****
 };
diff --git a/geom/geom/inc/TGeoOpticalSurface.h b/geom/geom/inc/TGeoOpticalSurface.h
new file mode 100644
index 00000000000..a56e6dc9218
--- /dev/null
+++ b/geom/geom/inc/TGeoOpticalSurface.h
@@ -0,0 +1,207 @@
+// @(#)root/geom:$Id$
+// Author: Andrei Gheata   05/12/18
+
+/*************************************************************************
+ * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers.               *
+ * All rights reserved.                                                  *
+ *                                                                       *
+ * For the licensing terms see $ROOTSYS/LICENSE.                         *
+ * For the list of contributors see $ROOTSYS/README/CREDITS.             *
+ *************************************************************************/
+
+#ifndef ROOT_TGeoOpticalSurface
+#define ROOT_TGeoOpticalSurface
+
+#include <TNamed.h>
+#include <TList.h>
+
+////////////////////////////////////////////////////////////////////////////
+//                                                                        //
+// TGeoOpticalSurface - class describing surface properties for           //
+//                      compatibility with Geant4                         //
+//                                                                        //
+////////////////////////////////////////////////////////////////////////////
+
+class TGeoOpticalSurface : public TNamed {
+public:
+   enum ESurfaceFinish {
+      kFpolished,             // smooth perfectly polished surface
+      kFpolishedfrontpainted, // smooth top-layer (front) paint
+      kFpolishedbackpainted,  // same is 'polished' but with a back-paint
+
+      kFground,             // rough surface
+      kFgroundfrontpainted, // rough top-layer (front) paint
+      kFgroundbackpainted,  // same as 'ground' but with a back-paint
+
+      kFpolishedlumirrorair,  // mechanically polished surface, with lumirror
+      kFpolishedlumirrorglue, // mechanically polished surface, with lumirror & meltmount
+      kFpolishedair,          // mechanically polished surface
+      kFpolishedteflonair,    // mechanically polished surface, with teflon
+      kFpolishedtioair,       // mechanically polished surface, with tio paint
+      kFpolishedtyvekair,     // mechanically polished surface, with tyvek
+      kFpolishedvm2000air,    // mechanically polished surface, with esr film
+      kFpolishedvm2000glue,   // mechanically polished surface, with esr film & meltmount
+
+      kFetchedlumirrorair,  // chemically etched surface, with lumirror
+      kFetchedlumirrorglue, // chemically etched surface, with lumirror & meltmount
+      kFetchedair,          // chemically etched surface
+      kFetchedteflonair,    // chemically etched surface, with teflon
+      kFetchedtioair,       // chemically etched surface, with tio paint
+      kFetchedtyvekair,     // chemically etched surface, with tyvek
+      kFetchedvm2000air,    // chemically etched surface, with esr film
+      kFetchedvm2000glue,   // chemically etched surface, with esr film & meltmount
+
+      kFgroundlumirrorair,  // rough-cut surface, with lumirror
+      kFgroundlumirrorglue, // rough-cut surface, with lumirror & meltmount
+      kFgroundair,          // rough-cut surface
+      kFgroundteflonair,    // rough-cut surface, with teflon
+      kFgroundtioair,       // rough-cut surface, with tio paint
+      kFgroundtyvekair,     // rough-cut surface, with tyvek
+      kFgroundvm2000air,    // rough-cut surface, with esr film
+      kFgroundvm2000glue,   // rough-cut surface, with esr film & meltmount
+
+      // for DAVIS model
+      kFRough_LUT,             // rough surface
+      kFRoughTeflon_LUT,       // rough surface wrapped in Teflon tape
+      kFRoughESR_LUT,          // rough surface wrapped with ESR
+      kFRoughESRGrease_LUT,    // rough surface wrapped with ESR and coupled with opical grease
+      kFPolished_LUT,          // polished surface
+      kFPolishedTeflon_LUT,    // polished surface wrapped in Teflon tape
+      kFPolishedESR_LUT,       // polished surface wrapped with ESR
+      kFPolishedESRGrease_LUT, // polished surface wrapped with ESR and coupled with opical grease
+      kFDetector_LUT           // polished surface with optical grease
+   };
+
+   enum ESurfaceModel {
+      kMglisur,  // original GEANT3 model
+      kMunified, // UNIFIED model
+      kMLUT,     // Look-Up-Table model
+      kMDAVIS,   // DAVIS model
+      kMdichroic // dichroic filter
+   };
+
+   enum ESurfaceType {
+      kTdielectric_metal,      // dielectric-metal interface
+      kTdielectric_dielectric, // dielectric-dielectric interface
+      kTdielectric_LUT,        // dielectric-Look-Up-Table interface
+      kTdielectric_LUTDAVIS,   // dielectric-Look-Up-Table DAVIS interface
+      kTdielectric_dichroic,   // dichroic filter interface
+      kTfirsov,                // for Firsov Process
+      kTx_ray                  // for x-ray mirror process
+   };
+
+private:
+   std::string fName = "";                  // Surface name
+   ESurfaceType fType = kTdielectric_metal; // Surface type
+   ESurfaceModel fModel = kMglisur;         // Surface model
+   ESurfaceFinish fFinish = kFpolished;     // Surface finish
+
+   Double_t fValue = 0.0;      // The value used to determine sigmaalpha and polish
+   Double_t fSigmaAlpha = 0.0; // The sigma of micro-facet polar angle
+   Double_t fPolish = 0.0;     // Polish parameter in glisur model
+
+   TList fProperties; // List of surface properties
+
+   // No copy
+   TGeoOpticalSurface(const TGeoOpticalSurface &);
+   TGeoOpticalSurface &operator=(const TGeoOpticalSurface &);
+
+public:
+   // constructors
+   TGeoOpticalSurface() {}
+
+   TGeoOpticalSurface(const char *name, ESurfaceModel model = kMglisur, ESurfaceFinish finish = kFpolished,
+                      ESurfaceType type = kTdielectric_dielectric, Double_t value = 1.0);
+
+   virtual ~TGeoOpticalSurface() {}
+
+   // Accessors
+   bool AddProperty(const char *property, const char *ref);
+   const char *GetPropertyRef(const char *property);
+   TList const &GetProperties() const { return fProperties; }
+   ESurfaceType GetType() const { return fType; }
+   ESurfaceModel GetModel() const { return fModel; }
+   ESurfaceFinish GetFinish() const { return fFinish; }
+   Double_t GetValue() const { return fValue; }
+   Double_t GetSigmaAlpha() const { return fSigmaAlpha; }
+
+   void SetType(ESurfaceType type) { fType = type; }
+   void SetModel(ESurfaceModel model) { fModel = model; }
+   void SetFinish(ESurfaceFinish finish) { fFinish = finish; }
+   void SetValue(Double_t value) { fValue = value; }
+   void SetSigmaAlpha(Double_t sigmaalpha) { fSigmaAlpha = sigmaalpha; }
+
+   void Print(Option_t *option = "") const;
+
+   static ESurfaceType StringToType(const char *type);
+   static const char *TypeToString(ESurfaceType type);
+   static ESurfaceModel StringToModel(const char *model);
+   static const char *ModelToString(ESurfaceModel model);
+   static ESurfaceFinish StringToFinish(const char *finish);
+   static const char *FinishToString(ESurfaceFinish finish);
+
+   ClassDef(TGeoOpticalSurface, 1) // Class representing an optical surface
+};
+
+////////////////////////////////////////////////////////////////////////////
+//                                                                        //
+// TGeoSkinSurface - class describing a surface having optical properties //
+//                      surrounding a volume                              //
+//                                                                        //
+////////////////////////////////////////////////////////////////////////////
+
+class TGeoVolume;
+
+class TGeoSkinSurface : public TNamed {
+private:
+   TGeoOpticalSurface const *fSurface = nullptr; // Referenced optical surface
+   TGeoVolume const *fVolume = nullptr;          // Referenced volume
+public:
+   TGeoSkinSurface() {}
+   TGeoSkinSurface(const char *name, const char *ref, TGeoOpticalSurface const *surf, TGeoVolume const *vol)
+      : TNamed(name, ref), fSurface(surf), fVolume(vol)
+   {
+   }
+   virtual ~TGeoSkinSurface() {}
+
+   TGeoOpticalSurface const *GetSurface() const { return fSurface; }
+   TGeoVolume const *GetVolume() const { return fVolume; }
+
+   void Print(Option_t *option = "") const;
+
+   ClassDef(TGeoSkinSurface, 1) // A surface with optical properties surrounding a volume
+};
+
+////////////////////////////////////////////////////////////////////////////
+//                                                                        //
+// TGeoBorderSurface - class describing a surface having optical          //
+//                      properties between 2 touching volumes             //
+//                                                                        //
+////////////////////////////////////////////////////////////////////////////
+
+class TGeoNode;
+
+class TGeoBorderSurface : public TNamed {
+private:
+   TGeoOpticalSurface const *fSurface = nullptr; // Referenced optical surface
+   TGeoNode const *fNode1 = nullptr;             // Referenced node 1
+   TGeoNode const *fNode2 = nullptr;             // Referenced node 2
+public:
+   TGeoBorderSurface() {}
+   TGeoBorderSurface(const char *name, const char *ref, TGeoOpticalSurface const *surf, TGeoNode const *node1,
+                     TGeoNode const *node2)
+      : TNamed(name, ref), fSurface(surf), fNode1(node1), fNode2(node2)
+   {
+   }
+   virtual ~TGeoBorderSurface() {}
+
+   TGeoOpticalSurface const *GetSurface() const { return fSurface; }
+   TGeoNode const *GetNode1() const { return fNode1; }
+   TGeoNode const *GetNode2() const { return fNode2; }
+
+   void Print(Option_t *option = "") const;
+
+   ClassDef(TGeoBorderSurface, 1) // A surface with optical properties betwqeen 2 touching volumes
+};
+
+#endif // ROOT_TGeoOpticalSurface
diff --git a/geom/geom/src/TGDMLMatrix.cxx b/geom/geom/src/TGDMLMatrix.cxx
new file mode 100644
index 00000000000..1f623daf068
--- /dev/null
+++ b/geom/geom/src/TGDMLMatrix.cxx
@@ -0,0 +1,99 @@
+// @(#)root/gdml:$Id$
+// Author: Andrei Gheata 05/12/2018
+
+/*************************************************************************
+ * Copyright (C) 1995-2011, Rene Brun and Fons Rademakers.               *
+ * All rights reserved.                                                  *
+ *                                                                       *
+ * For the licensing terms see $ROOTSYS/LICENSE.                         *
+ * For the list of contributors see $ROOTSYS/README/CREDITS.             *
+ *************************************************************************/
+
+/** \class TGDMLMatrix
+\ingroup Geometry_gdml
+  This class is used in the process of reading and writing the GDML "matrix" tag.
+It represents a matrix with arbitrary number of rows and columns, storing elements
+in double precision.
+*/
+
+#include "TGDMLMatrix.h"
+
+#include <cassert>
+
+ClassImp(TGDMLMatrix)
+
+//_____________________________________________________________________________
+TGDMLMatrix::TGDMLMatrix(const char *name, size_t rows,size_t cols)
+   : TNamed(name, "")
+{
+// Constructor
+   if ((rows <= 0) || (cols <= 0))
+   {
+     Fatal("TGDMLMatrix::TGDMLMatrix(rows,cols)", "Wrong number of rows/cols");
+   }
+   fNrows = rows;
+   fNcols = cols;
+   fNelem = rows * cols;
+   fMatrix = new Double_t[fNelem];
+}
+
+//_____________________________________________________________________________
+TGDMLMatrix::TGDMLMatrix(const TGDMLMatrix& rhs)
+  : TNamed(rhs), fNelem(rhs.fNelem), fNrows(rhs.fNrows), fNcols(rhs.fNcols), fMatrix(nullptr)
+{
+// Copy constructor
+   if (rhs.fMatrix)
+   {
+     fMatrix = new Double_t[fNelem];
+     memcpy(fMatrix, rhs.fMatrix, fNelem * sizeof(Double_t));
+   }
+}
+
+//_____________________________________________________________________________
+TGDMLMatrix& TGDMLMatrix::operator=(const TGDMLMatrix& rhs)
+{
+// Assignment
+   if (this == &rhs)  { return *this; }
+   TNamed::operator=(rhs);
+   fNrows = rhs.fNrows;
+   fNcols = rhs.fNcols;
+   fNelem = fNrows * fNcols;
+   if (rhs.fMatrix)
+   {
+     fMatrix = new Double_t[fNelem];
+     memcpy(fMatrix, rhs.fMatrix, fNelem * sizeof(Double_t));
+   }
+   return *this;
+}
+
+//_____________________________________________________________________________
+void TGDMLMatrix::Set(size_t r, size_t c, Double_t a)
+{ 
+   assert(r < fNrows && c < fNcols);
+   fMatrix[fNcols*r+c] = a;
+}
+
+//_____________________________________________________________________________
+Double_t TGDMLMatrix::Get(size_t r, size_t c) const
+{
+   assert(r < fNrows && c < fNcols);
+   return fMatrix[fNcols*r+c];
+}
+
+//_____________________________________________________________________________
+void TGDMLMatrix::Print(Option_t *) const
+{
+// Print info about this matrix
+   printf("*** matrix: %-20s coldim = %ld  rows = %ld\n", GetName(), fNcols, fNrows);
+   if (!fTitle.IsNull()) {
+      printf("   %s\n", fTitle.Data());
+      return;
+   }
+   for (size_t row = 0; row < fNrows; ++row) {
+         printf("   ");
+      for (size_t col = 0; col < fNcols; ++col) {
+         printf("%8.3g", Get(row, col));
+      }
+      printf("\n");
+   }
+}
diff --git a/geom/geom/src/TGeoManager.cxx b/geom/geom/src/TGeoManager.cxx
index a656d60a8e6..121f6832c53 100644
--- a/geom/geom/src/TGeoManager.cxx
+++ b/geom/geom/src/TGeoManager.cxx
@@ -283,6 +283,8 @@ in order to enhance rays.
 #include "TEnv.h"
 #include "TGeoParallelWorld.h"
 #include "TGeoRegion.h"
+#include "TGDMLMatrix.h"
+#include "TGeoOpticalSurface.h"
 
 // statics and globals
 
@@ -341,7 +343,12 @@ TGeoManager::TGeoManager()
       fNtracks = 0;
       fNpdg = 0;
       fPdgNames = 0;
+      fGDMLMatrices = 0;
+      fOpticalSurfaces = 0;
+      fSkinSurfaces = 0;
+      fBorderSurfaces = 0;
       memset(fPdgId, 0, 1024*sizeof(Int_t));
+//   TObjArray            *fNavigators;       //! list of navigators
       fCurrentTrack = 0;
       fCurrentVolume = 0;
       fTopVolume = 0;
@@ -444,6 +451,10 @@ void TGeoManager::Init()
    fNtracks = 0;
    fNpdg = 0;
    fPdgNames = 0;
+   fGDMLMatrices = new TObjArray();
+   fOpticalSurfaces = new TObjArray();
+   fSkinSurfaces = new TObjArray();
+   fBorderSurfaces = new TObjArray();
    memset(fPdgId, 0, 1024*sizeof(Int_t));
    fCurrentTrack = 0;
    fCurrentVolume = 0;
@@ -523,6 +534,10 @@ TGeoManager::TGeoManager(const TGeoManager& gm) :
   fGVolumes(gm.fGVolumes),
   fTracks(gm.fTracks),
   fPdgNames(gm.fPdgNames),
+  fGDMLMatrices(gm.fGDMLMatrices),
+  fOpticalSurfaces(gm.fOpticalSurfaces),
+  fSkinSurfaces(gm.fSkinSurfaces),
+  fBorderSurfaces(gm.fBorderSurfaces),
   fMaterials(gm.fMaterials),
   fMedia(gm.fMedia),
   fNodes(gm.fNodes),
@@ -608,6 +623,10 @@ TGeoManager& TGeoManager::operator=(const TGeoManager& gm)
       fGVolumes=gm.fGVolumes;
       fTracks=gm.fTracks;
       fPdgNames=gm.fPdgNames;
+      fGDMLMatrices=gm.fGDMLMatrices;
+      fOpticalSurfaces = gm.fOpticalSurfaces;
+      fSkinSurfaces = gm.fSkinSurfaces;
+      fBorderSurfaces = gm.fBorderSurfaces;
       fMaterials=gm.fMaterials;
       fMedia=gm.fMedia;
       fNodes=gm.fNodes;
@@ -684,6 +703,10 @@ TGeoManager::~TGeoManager()
    if (fTracks) {fTracks->Delete(); SafeDelete( fTracks );}
    SafeDelete( fUniqueVolumes );
    if (fPdgNames) {fPdgNames->Delete(); SafeDelete( fPdgNames );}
+   if (fGDMLMatrices) {fGDMLMatrices->Delete(); SafeDelete( fGDMLMatrices );}
+   if (fOpticalSurfaces) {fOpticalSurfaces->Delete(); SafeDelete( fOpticalSurfaces );}
+   if (fSkinSurfaces) {fSkinSurfaces->Delete(); SafeDelete( fSkinSurfaces );}
+   if (fBorderSurfaces) {fBorderSurfaces->Delete(); SafeDelete( fBorderSurfaces );}
    ClearNavigators();
    CleanGarbage();
    SafeDelete( fPainter );
@@ -1928,6 +1951,82 @@ void TGeoManager::SetPdgName(Int_t pdg, const char *name)
    fPdgNames->AddAtAndExpand(pdgname, fNpdg++);
 }
 
+////////////////////////////////////////////////////////////////////////////////
+/// Get GDML matrix with a given name;
+
+TGDMLMatrix *TGeoManager::GetGDMLMatrix(const char *name) const
+{
+   return (TGDMLMatrix*)fGDMLMatrices->FindObject(name);
+}
+
+////////////////////////////////////////////////////////////////////////////////
+/// Add GDML matrix;
+void TGeoManager::AddGDMLMatrix(TGDMLMatrix *mat)
+{
+   if (GetGDMLMatrix(mat->GetName())) {
+      Error("AddGDMLMatrix", "Matrix %s already added to manager", mat->GetName());
+      return;
+   }
+   fGDMLMatrices->Add(mat);
+}
+
+////////////////////////////////////////////////////////////////////////////////
+/// Get optical surface with a given name;
+
+TGeoOpticalSurface *TGeoManager::GetOpticalSurface(const char *name) const
+{
+   return (TGeoOpticalSurface*)fOpticalSurfaces->FindObject(name);
+}
+
+////////////////////////////////////////////////////////////////////////////////
+/// Add optical surface;
+void TGeoManager::AddOpticalSurface(TGeoOpticalSurface *optsurf)
+{
+   if (GetOpticalSurface(optsurf->GetName())) {
+      Error("AddOpticalSurface", "Surface %s already added to manager", optsurf->GetName());
+      return;
+   }
+   fOpticalSurfaces->Add(optsurf);
+}
+
+////////////////////////////////////////////////////////////////////////////////
+/// Get skin surface with a given name;
+
+TGeoSkinSurface *TGeoManager::GetSkinSurface(const char *name) const
+{
+   return (TGeoSkinSurface*)fSkinSurfaces->FindObject(name);
+}
+
+////////////////////////////////////////////////////////////////////////////////
+/// Add skin surface;
+void TGeoManager::AddSkinSurface(TGeoSkinSurface *surf)
+{
+   if (GetSkinSurface(surf->GetName())) {
+      Error("AddSkinSurface", "Surface %s already added to manager", surf->GetName());
+      return;
+   }
+   fSkinSurfaces->Add(surf);
+}
+
+////////////////////////////////////////////////////////////////////////////////
+/// Get border surface with a given name;
+
+TGeoBorderSurface *TGeoManager::GetBorderSurface(const char *name) const
+{
+   return (TGeoBorderSurface*)fBorderSurfaces->FindObject(name);
+}
+
+////////////////////////////////////////////////////////////////////////////////
+/// Add border surface;
+void TGeoManager::AddBorderSurface(TGeoBorderSurface *surf)
+{
+   if (GetBorderSurface(surf->GetName())) {
+      Error("AddBorderSurface", "Surface %s already added to manager", surf->GetName());
+      return;
+   }
+   fBorderSurfaces->Add(surf);
+}
+
 ////////////////////////////////////////////////////////////////////////////////
 /// Fill node copy numbers of current branch into an array.
 
diff --git a/geom/geom/src/TGeoMaterial.cxx b/geom/geom/src/TGeoMaterial.cxx
index 9bcd1de35cb..c3207756d29 100644
--- a/geom/geom/src/TGeoMaterial.cxx
+++ b/geom/geom/src/TGeoMaterial.cxx
@@ -225,6 +225,10 @@ TGeoMaterial::TGeoMaterial(const TGeoMaterial& gm) :
 
 {
    //copy constructor
+   fProperties.SetOwner();
+   TIter next(&fProperties);
+   TNamed *property;
+   while ((property = (TNamed*)next())) fProperties.Add(new TNamed(*property));
 }
 
 ////////////////////////////////////////////////////////////////////////////////
@@ -249,6 +253,10 @@ TGeoMaterial& TGeoMaterial::operator=(const TGeoMaterial& gm)
       fElement=gm.fElement;
       fUserExtension = gm.fUserExtension->Grab();
       fFWExtension = gm.fFWExtension->Grab();
+      fProperties.SetOwner();
+      TIter next(&fProperties);
+      TNamed *property;
+      while ((property = (TNamed*)next())) fProperties.Add(new TNamed(*property));
    }
    return *this;
 }
@@ -277,6 +285,27 @@ void TGeoMaterial::SetUserExtension(TGeoExtension *ext)
    if (ext) fUserExtension = ext->Grab();
 }
 
+//_____________________________________________________________________________
+const char *TGeoMaterial::GetPropertyRef(const char *property)
+{
+   // Find reference for a given property
+   TNamed *prop = (TNamed*)fProperties.FindObject(property);
+   return (prop) ? prop->GetTitle() : nullptr;
+}
+
+//_____________________________________________________________________________
+bool TGeoMaterial::AddProperty(const char *property, const char *ref)
+{
+   fProperties.SetOwner();
+   if (GetPropertyRef(property)) {
+      Error("AddProperty", "Property %s already added to material %s",
+         property, GetName());
+      return false;
+   }
+   fProperties.Add(new TNamed(property, ref));
+   return true;
+}
+
 ////////////////////////////////////////////////////////////////////////////////
 /// Connect framework defined extension to the material. The material "grabs" a copy,
 /// so the original object can be released by the producer. Release the previously
diff --git a/geom/geom/src/TGeoOpticalSurface.cxx b/geom/geom/src/TGeoOpticalSurface.cxx
new file mode 100644
index 00000000000..3ba9ca6ee12
--- /dev/null
+++ b/geom/geom/src/TGeoOpticalSurface.cxx
@@ -0,0 +1,300 @@
+// @(#)root/geom:$Id$
+// Author: Andrei Gheata   05/12/18
+
+/*************************************************************************
+ * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers.               *
+ * All rights reserved.                                                  *
+ *                                                                       *
+ * For the licensing terms see $ROOTSYS/LICENSE.                         *
+ * For the list of contributors see $ROOTSYS/README/CREDITS.             *
+ *************************************************************************/
+
+
+/** \class TGeoMedium
+\ingroup Geometry_classes
+
+This is a wrapper class to G4OpticalSurface
+*/
+
+#include "TGeoOpticalSurface.h"
+
+#include <string>
+
+#include "TGeoVolume.h"
+#include "TGeoNode.h"
+
+ClassImp(TGeoOpticalSurface)
+ClassImp(TGeoSkinSurface)
+ClassImp(TGeoBorderSurface)
+
+//_____________________________________________________________________________
+TGeoOpticalSurface::TGeoOpticalSurface(const char *name, ESurfaceModel model, ESurfaceFinish finish,
+                                       ESurfaceType type, Double_t value)
+   : TNamed(name, ""), fType(type), fModel(model), fFinish(finish), fValue(value)
+{
+   // Constructor
+   fProperties.SetOwner();
+   if (model == kMglisur) {
+      fPolish = value;
+      fSigmaAlpha = 0.0;
+   } else if (model == kMunified || model == kMLUT || model == kMdichroic || model == kMDAVIS) {
+      fSigmaAlpha = value;
+      fPolish = 0.0;
+   } else {
+      Fatal("TGeoOpticalSurface::TGeoOpticalSurface()", "Constructor called with INVALID model.");
+   }
+}
+
+//_____________________________________________________________________________
+TGeoOpticalSurface::ESurfaceType TGeoOpticalSurface::StringToType(const char *name)
+{
+   // Convert string to optical surface type
+   ESurfaceType type;
+   TString stype(name);
+   if ((stype == "dielectric_metal") || (stype == "0")) {
+      type = kTdielectric_metal;
+   } else if ((stype == "dielectric_dielectric") || (stype == "1")) {
+      type = kTdielectric_dielectric;
+   } else if ((stype == "dielectric_LUT") || (stype == "2")) {
+      type = kTdielectric_LUT;
+   } else if ((stype == "dielectric_dichroic") || (stype == "3")) {
+      type = kTdielectric_dichroic;
+   } else if ((stype == "firsov") || (stype == "4")) {
+      type = kTfirsov;
+   } else {
+      type = kTx_ray;
+   }
+   return type;
+}
+
+//_____________________________________________________________________________
+const char *TGeoOpticalSurface::TypeToString(ESurfaceType type)
+{
+   // Convert surface type to string
+   switch (type) {
+   case kTdielectric_metal: return "dielectric_metal";
+   case kTdielectric_dielectric: return "dielectric_dielectric";
+   case kTdielectric_LUT: return "dielectric_LUT";
+   case kTdielectric_dichroic: return "dielectric_dichroic";
+   case kTfirsov: return "firsov";
+   case kTx_ray: return "x_ray";
+   case kTdielectric_LUTDAVIS:;
+   }
+
+   return "unhandled surface type";
+}
+
+//_____________________________________________________________________________
+TGeoOpticalSurface::ESurfaceModel TGeoOpticalSurface::StringToModel(const char *name)
+{
+
+   // Convert string to optical surface type
+   TString smodel(name);
+   ESurfaceModel model;
+   if ((smodel == "glisur") || (smodel == "0")) {
+      model = kMglisur;
+   } else if ((smodel == "unified") || (smodel == "1")) {
+      model = kMunified;
+   } else if ((smodel == "LUT") || (smodel == "2")) {
+      model = kMLUT;
+   } else {
+      model = kMdichroic;
+   }
+   return model;
+}
+
+//_____________________________________________________________________________
+const char *TGeoOpticalSurface::ModelToString(ESurfaceModel model)
+{
+   // Convert optical surface model to string
+   switch (model) {
+   case kMglisur: return "glisur";
+   case kMunified: return "unified";
+   case kMLUT: return "LUT";
+   case kMdichroic: return "dichoic";
+   case kMDAVIS:;
+   }
+
+   return "unhandled model type";
+}
+
+//_____________________________________________________________________________
+TGeoOpticalSurface::ESurfaceFinish TGeoOpticalSurface::StringToFinish(const char *name)
+{
+   // Convert surface finish to string
+   TString sfinish(name);
+   ESurfaceFinish finish;
+   if ((sfinish == "polished") || (sfinish == "0")) {
+      finish = kFpolished;
+   } else if ((sfinish == "polishedfrontpainted") || (sfinish == "1")) {
+      finish = kFpolishedfrontpainted;
+   } else if ((sfinish == "polishedbackpainted") || (sfinish == "2")) {
+      finish = kFpolishedbackpainted;
+   } else if ((sfinish == "ground") || (sfinish == "3")) {
+      finish = kFground;
+   } else if ((sfinish == "groundfrontpainted") || (sfinish == "4")) {
+      finish = kFgroundfrontpainted;
+   } else if ((sfinish == "groundbackpainted") || (sfinish == "5")) {
+      finish = kFgroundbackpainted;
+   } else if ((sfinish == "polishedlumirrorair") || (sfinish == "6")) {
+      finish = kFpolishedlumirrorair;
+   } else if ((sfinish == "polishedlumirrorglue") || (sfinish == "7")) {
+      finish = kFpolishedlumirrorglue;
+   } else if ((sfinish == "polishedair") || (sfinish == "8")) {
+      finish = kFpolishedair;
+   } else if ((sfinish == "polishedteflonair") || (sfinish == "9")) {
+      finish = kFpolishedteflonair;
+   } else if ((sfinish == "polishedtioair") || (sfinish == "10")) {
+      finish = kFpolishedtioair;
+   } else if ((sfinish == "polishedtyvekair") || (sfinish == "11")) {
+      finish = kFpolishedtyvekair;
+   } else if ((sfinish == "polishedvm2000air") || (sfinish == "12")) {
+      finish = kFpolishedvm2000air;
+   } else if ((sfinish == "polishedvm2000glue") || (sfinish == "13")) {
+      finish = kFpolishedvm2000glue;
+   } else if ((sfinish == "etchedlumirrorair") || (sfinish == "14")) {
+      finish = kFetchedlumirrorair;
+   } else if ((sfinish == "etchedlumirrorglue") || (sfinish == "15")) {
+      finish = kFetchedlumirrorglue;
+   } else if ((sfinish == "etchedair") || (sfinish == "16")) {
+      finish = kFetchedair;
+   } else if ((sfinish == "etchedteflonair") || (sfinish == "17")) {
+      finish = kFetchedteflonair;
+   } else if ((sfinish == "etchedtioair") || (sfinish == "18")) {
+      finish = kFetchedtioair;
+   } else if ((sfinish == "etchedtyvekair") || (sfinish == "19")) {
+      finish = kFetchedtyvekair;
+   } else if ((sfinish == "etchedvm2000air") || (sfinish == "20")) {
+      finish = kFetchedvm2000air;
+   } else if ((sfinish == "etchedvm2000glue") || (sfinish == "21")) {
+      finish = kFetchedvm2000glue;
+   } else if ((sfinish == "groundlumirrorair") || (sfinish == "22")) {
+      finish = kFgroundlumirrorair;
+   } else if ((sfinish == "groundlumirrorglue") || (sfinish == "23")) {
+      finish = kFgroundlumirrorglue;
+   } else if ((sfinish == "groundair") || (sfinish == "24")) {
+      finish = kFgroundair;
+   } else if ((sfinish == "groundteflonair") || (sfinish == "25")) {
+      finish = kFgroundteflonair;
+   } else if ((sfinish == "groundtioair") || (sfinish == "26")) {
+      finish = kFgroundtioair;
+   } else if ((sfinish == "groundtyvekair") || (sfinish == "27")) {
+      finish = kFgroundtyvekair;
+   } else if ((sfinish == "groundvm2000air") || (sfinish == "28")) {
+      finish = kFgroundvm2000air;
+   } else {
+      finish = kFgroundvm2000glue;
+   }
+
+   return finish;
+}
+
+//_____________________________________________________________________________
+const char *TGeoOpticalSurface::FinishToString(ESurfaceFinish finish)
+{
+   switch (finish) {
+   case kFpolished: return "polished";
+   case kFpolishedfrontpainted: return "polishedfrontpainted";
+   case kFpolishedbackpainted: return "polishedbackpainted";
+
+   case kFground: return "ground";
+   case kFgroundfrontpainted: return "groundfrontpainted";
+   case kFgroundbackpainted: return "groundbackpainted";
+
+   case kFpolishedlumirrorair: return "polishedlumirrorair";
+   case kFpolishedlumirrorglue: return "polishedlumirrorglue";
+   case kFpolishedair: return "polishedair";
+   case kFpolishedteflonair: return "polishedteflonair";
+   case kFpolishedtioair: return "polishedtioair";
+   case kFpolishedtyvekair: return "polishedtyvekair";
+   case kFpolishedvm2000air: return "polishedvm2000air";
+   case kFpolishedvm2000glue: return "polishedvm2000glue";
+
+   case kFetchedlumirrorair: return "etchedlumirrorair";
+   case kFetchedlumirrorglue: return "etchedlumirrorglue";
+   case kFetchedair: return "etchedair";
+   case kFetchedteflonair: return "etchedteflonair";
+   case kFetchedtioair: return "etchedtioair";
+   case kFetchedtyvekair: return "etchedtyvekair";
+   case kFetchedvm2000air: return "etchedvm2000air";
+   case kFetchedvm2000glue: return "etchedvm2000glue";
+
+   case kFgroundlumirrorair: return "groundlumirrorair";
+   case kFgroundlumirrorglue: return "groundlumirrorglue";
+   case kFgroundair: return "groundair";
+   case kFgroundteflonair: return "groundteflonair";
+   case kFgroundtioair: return "groundtioair";
+   case kFgroundtyvekair: return "groundtyvekair";
+   case kFgroundvm2000air: return "groundvm2000air";
+   case kFgroundvm2000glue: return "groundvm2000glue";
+   case kFRough_LUT:
+   case kFRoughTeflon_LUT:
+   case kFRoughESR_LUT:
+   case kFRoughESRGrease_LUT:
+   case kFPolished_LUT:
+   case kFPolishedTeflon_LUT:
+   case kFPolishedESR_LUT:
+   case kFPolishedESRGrease_LUT:
+   case kFDetector_LUT:;
+   }
+
+   return "unhandled model finish";
+}
+
+//_____________________________________________________________________________
+const char *TGeoOpticalSurface::GetPropertyRef(const char *property)
+{
+   // Find reference for a given property
+   TNamed *prop = (TNamed *)fProperties.FindObject(property);
+   return (prop) ? prop->GetTitle() : nullptr;
+}
+
+//_____________________________________________________________________________
+bool TGeoOpticalSurface::AddProperty(const char *property, const char *ref)
+{
+   fProperties.SetOwner();
+   if (GetPropertyRef(property)) {
+      Error("AddProperty", "Property %s already added to optical surface %s", property, GetName());
+      return false;
+   }
+   fProperties.Add(new TNamed(property, ref));
+   return true;
+}
+
+//_____________________________________________________________________________
+void TGeoOpticalSurface::Print(Option_t *) const
+{
+   // Print info about this optical surface
+   printf("*** opticalsurface: %s   type: %s   model: %s   finish: %s   value = %g\n", GetName(),
+          TGeoOpticalSurface::TypeToString(fType), TGeoOpticalSurface::ModelToString(fModel),
+          TGeoOpticalSurface::FinishToString(fFinish), fValue);
+   if (fProperties.GetSize()) {
+      TIter next(&fProperties);
+      TNamed *property;
+      while ((property = (TNamed *)next()))
+         printf("   property: %s ref: %s\n", property->GetName(), property->GetTitle());
+   }
+}
+
+//_____________________________________________________________________________
+void TGeoSkinSurface::Print(Option_t *) const
+{
+   // Print info about this optical surface
+   if (!fVolume) {
+      Error("Print", "Skin surface %s: volume not set", GetName());
+      return;
+   }
+   printf("*** skinsurface: %s   surfaceproperty: %s   volumeref: %s \n", GetName(), GetTitle(), fVolume->GetName());
+}
+
+//_____________________________________________________________________________
+void TGeoBorderSurface::Print(Option_t *) const
+{
+   // Print info about this optical surface
+   if (!fNode1 || !fNode2) {
+      Error("Print", "Border surface %s: nodes not set", GetName());
+      return;
+   }
+   printf("*** bordersurface: %s   surfaceproperty: %s   physvolref: %s  %s \n", GetName(), GetTitle(),
+          fNode1->GetName(), fNode2->GetName());
+}
diff --git a/tutorials/geom/gdml/opticalsurfaces.gdml b/tutorials/geom/gdml/opticalsurfaces.gdml
new file mode 100644
index 00000000000..8d4fd72fb5c
--- /dev/null
+++ b/tutorials/geom/gdml/opticalsurfaces.gdml
@@ -0,0 +1,78 @@
+<?xml version="1.0" encoding="UTF-8" ?>
+<gdml xmlns:gdml="http://cern.ch/2001/Schemas/GDML" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:noNamespaceSchemaLocation="http://service-spi.web.cern.ch/service-spi/app/releases/GDML/schema/gdml.xsd" >
+
+ <define>
+  <position name="detinWorldp" unit="mm" x="0" y="0" z="0" />
+  <position name="pos2" unit="mm" x="2000" y="0" z="0" /> 
+
+  <matrix name="prop1" coldim="2" values="1.0 7.0 
+                                          2.0 9.4"/>
+
+  <matrix name="prop2" coldim="2" values="3.0 -1.0 
+                                          4.0 -9.4
+                                          5.0 -11.1"/>
+
+  <matrix name="RINDEX" coldim="2" values="1.65*eV 1.58" />
+  <matrix name="REFLECTIVITY" coldim="2" values="1.65*eV 0.8" />
+ </define>
+
+ <materials>
+  <element Z="7" formula="N" name="Nitrogen" >
+   <atom value="14.01" />
+  </element>
+  <element Z="8" formula="O" name="Oxygen" >
+   <atom value="16" />
+  </element>
+
+  <material formula=" " name="Air" >
+   <property name="myproperty1" ref="prop1"/>
+   <property name="myproperty2" ref="prop2"/>
+   <property name="RINDEX" ref="RINDEX"/>  
+   <D value="0.00129" />
+   <fraction n="0.7" ref="Nitrogen" />
+   <fraction n="0.3" ref="Oxygen" />
+  </material>
+ </materials>
+
+ <solids>
+  <box aunit="radian" lunit="mm" name="world" x="10000" y="10000" z="10000" />
+  <box aunit="radian" lunit="mm" name="det" x="2000" y="2000" z="2000" />
+  <opticalsurface name="surf1" model="glisur" finish="polished" type="dielectric_dielectric" value="1.0"/>
+  <opticalsurface name="surf2" model="glisur" finish="polished" type="dielectric_metal" value="1.0">
+    <property name="REFLECTIVITY" ref="REFLECTIVITY" />
+  </opticalsurface>
+ </solids>
+
+ <structure>
+  <volume name="Detector" >
+   <materialref ref="Air" />
+   <solidref ref="det" />
+  </volume>
+  <volume name="World" >
+   <materialref ref="Air" />
+   <solidref ref="world" />
+   <physvol name="pv1">
+    <volumeref ref="Detector" />
+    <positionref ref="detinWorldp" />
+   </physvol>
+   <physvol name="pv2">
+    <volumeref ref="Detector" />
+    <positionref ref="pos2" />
+   </physvol>
+  </volume>
+
+  <skinsurface name="skinsrf1" surfaceproperty="surf1" >
+    <volumeref ref="Detector"/>
+  </skinsurface> 
+
+  <bordersurface name="bordersrf1" surfaceproperty="surf2" >
+    <physvolref ref="pv1"/>
+    <physvolref ref="pv2"/>
+  </bordersurface> 
+
+ </structure>
+
+ <setup name="Default" version="1.0" >
+  <world ref="World" />
+ </setup>
+</gdml>
diff --git a/tutorials/geom/gdml/testoptical.C b/tutorials/geom/gdml/testoptical.C
new file mode 100644
index 00000000000..273fedb0e95
--- /dev/null
+++ b/tutorials/geom/gdml/testoptical.C
@@ -0,0 +1,69 @@
+/// \file
+/// \ingroup tutorial_geom
+/// Tests importing/exporting optical surfaces from GDML
+///
+/// Optical surfaces, skin surfaces and border surfaces are imported in object arrays
+/// stored by TGeoManager class. Optical surfaces do not store property arrays but point
+/// to GDML matrices describing such properties. One can get the data for such property
+/// like:
+///   TGeoOpticalSurface *surf = geom->GetOpticalSurface("surf1");
+///   const char *property = surf=>GetPropertyRef("REFLECTIVITY");
+///   TGeoGDMLMatrix *m = geom->GetGDMLMatrix(property);
+/// Skin surfaces and border surfaces can be retrieved from the TGeoManager object by using:
+///   TObjArray *skin_array = geom->GetListOfSkinSurfaces();
+///   TObjArra8 *border_array = geom->GetListOfBorderSurfaces();
+/// Alternatively accessors by name can also be used: GetSkinSurface(name)/GetBorderSurface(name)
+///
+/// \author Andrei Gheata
+
+#include <cassert>
+#include <TObjArray.h>
+#include <TROOT.h>
+#include <TGeoOpticalSurface.h>
+#include <TGeoManager.h>
+
+double Checksum(TGeoManager *geom)
+{
+   double sum = 0.;
+   TIter next(geom->GetListOfOpticalSurfaces());
+   TGeoOpticalSurface *surf;
+   while ((surf = (TGeoOpticalSurface *)next())) {
+      sum += (double)surf->GetType() + (double)surf->GetModel() + (double)surf->GetFinish() + surf->GetValue();
+      TString name = surf->GetName();
+      sum += (double)name.Hash();
+      name = surf->GetTitle();
+      sum += (double)name.Hash();
+   }
+   return sum;
+}
+
+int testoptical()
+{
+   TString geofile = gROOT->GetTutorialDir() + "/geom/gdml/opticalsurfaces.gdml";
+   TGeoManager::SetExportPrecision(8);
+   TGeoManager::SetVerboseLevel(0);
+   printf("=== Importing %s ...\n", geofile.Data());
+   TGeoManager *geom = TGeoManager::Import(geofile);
+   printf("=== List of GDML matrices:\n");
+   geom->GetListOfGDMLMatrices()->Print();
+   printf("=== List of optical surfaces:\n");
+   geom->GetListOfOpticalSurfaces()->Print();
+   printf("=== List of skin surfaces:\n");
+   geom->GetListOfSkinSurfaces()->Print();
+   printf("=== List of border surfaces:\n");
+   geom->GetListOfBorderSurfaces()->Print();
+   // Compute some checksum for optical surfaces
+   double checksum1 = Checksum(geom);
+   printf("=== Exporting as .gdml, then importing back\n");
+   geom->Export("tmp.gdml");
+   geom = TGeoManager::Import("tmp.gdml");
+   double checksum2 = Checksum(geom);
+   assert((checksum2 == checksum1) && "Exporting/importing as .gdml not OK");
+   printf("=== Exporting as .root, then importing back\n");
+   geom->Export("tmp.root");
+   geom = TGeoManager::Import("tmp.root");
+   double checksum3 = Checksum(geom);
+   assert((checksum3 == checksum1) && "Exporting/importing as .root not OK");
+   printf("all OK\n");
+   return 0;
+}
-- 
GitLab