Fix creation of tessellated shapes from CAD files.

DDCAD/include/DDCAD/Utilities.h

@@ -15,6 +15,7 @@
 
 #include <vector>
 #include <cmath>
+#include <limits>
 
 #include <TGeoTessellated.h>
 #include <TGeoVector3.h>
@@ -25,52 +26,66 @@ namespace dd4hep {
   /// Namespace for implementation details of the AIDA detector description toolkit
   namespace cad  {
 
-    inline std::stringstream streamFacet(TGeoFacet const& facet,
-					 TGeoTessellated const& shape) {
+#if ROOT_VERSION_CODE >= ROOT_VERSION(6,31,1)
+    inline std::string streamFacet(TGeoFacet const& facet,
+                                   TGeoTessellated const& shape) {
       using ::operator<<;
       std::stringstream str;
-#if ROOT_VERSION_CODE >= ROOT_VERSION(6,31,1)
       str << "{";
       for (int i = 0; i < facet.GetNvert(); ++i) {
-	str << shape.GetVertex(facet[i]);
-	if (i != facet.GetNvert() - 1)
-	  str << ", ";
+        str << shape.GetVertex(facet[i]);
+        if (i != facet.GetNvert() - 1)
+          str << ", ";
       }
       str << "}";
+      return str.str();
+    }
 #else
+    inline std::string streamFacet(TGeoFacet const& facet,
+                                   TGeoTessellated const& /* shape */) {
+      using ::operator<<;
+      std::stringstream str;
       str << facet;
-#endif
-      return str;
+      return str.str();
     }
-
-    inline std::stringstream streamVertices(ROOT::Geom::Vertex_t const& v1,
+#endif
+
+    inline std::string streamVertices(ROOT::Geom::Vertex_t const& v1,
 					    ROOT::Geom::Vertex_t const& v2,
 					    ROOT::Geom::Vertex_t const& v3) {
       using ::operator<<;
       std::stringstream str;
       str << "{" << v1 << ", " << v2 <<  ", " << v3 << "}";
-      return str;
+      return str.str();
     }
-
     // Determine if the facet is degenerated by calculating its determinant
     inline bool facetIsDegenerated(std::vector<ROOT::Geom::Vertex_t> const& vertices){
-      const ROOT::Geom::Vertex_t& v1 = vertices[0];
-      const ROOT::Geom::Vertex_t& v2 = vertices[1];
-      const ROOT::Geom::Vertex_t& v3 = vertices[2];
-      constexpr double epsilon = 1.e-20;
-      // v1.x v2.x v3.x v1.x v2.x
-      //
-      // v1.y v2.y v3.y	v1.y v2.y
-      //
-      // v1.z v2.z v3.z	v1.z v2.z
-      double det =  0.0
-        + v1.x() * v2.y() * v3.z()
-        + v2.x() * v3.y() * v1.z()
-        + v3.x() * v1.y() * v2.z()
-        - v1.z() * v2.y() * v3.x()
-        - v2.z() * v3.y() * v1.x()
-        - v3.z() * v1.y() * v2.x();
-      return std::fabs(det) < epsilon;
+      using ROOT::Geom::Vertex_t;
+      // Compute normal using non-zero segments
+      constexpr double kTolerance = 1.e-20;
+      bool degenerated = true;
+      Vertex_t normal;
+      int fNvert = int(vertices.size());
+      for (int i = 0; i < fNvert - 1; ++i) {
+        Vertex_t e1 = vertices[i + 1] - vertices[i];
+        if (e1.Mag2() < kTolerance)
+          continue;
+        for (int j = i + 1; j < fNvert; ++j) {
+          Vertex_t e2 = vertices[(j + 1) % fNvert] - vertices[j];
+          if (e2.Mag2() < kTolerance)
+            continue;
+          normal = Vertex_t::Cross(e1, e2);
+          // e1 and e2 may be colinear
+          if (normal.Mag2() < kTolerance)
+            continue;
+          normal.Normalize();
+          degenerated = false;
+          break;
+        }
+        if (!degenerated)
+          break;
+      }
+      return degenerated;
     }
   }
 }

DDCAD/src/ASSIMPReader.cpp

@@ -63,9 +63,9 @@ ASSIMPReader::readShapes(const std::string& source, double unit_length)  const
         if ( dump_facets )   {
           for( size_t i=0, n=shape->GetNfacets(); i < n; ++i )   {
             const auto& facet = shape->GetFacet(i);
-            std::stringstream str = dd4hep::cad::streamFacet(facet, shape);
+            std::string str = dd4hep::cad::streamFacet(facet, shape);
             printout(ALWAYS,"ASSIMPReader","++ Facet %4ld : %s",
-                     i, str.str().c_str());
+                     i, str.c_str());
           }
         }
         shape->SetTitle(TESSELLATEDSOLID_TAG);
@@ -90,7 +90,6 @@ ASSIMPReader::readVolumes(const std::string& source, double unit_length)  const
   bool dump_facets = ((flags>>8)&0x1) == 1;
   int aiflags = aiProcess_Triangulate|aiProcess_JoinIdenticalVertices|aiProcess_CalcTangentSpace;
   auto scene = importer->ReadFile( source.c_str(), aiflags);
-  char text[1048];
 
   if ( !scene )  {
     except("ASSIMPReader","+++ FileNotFound: %s",source.c_str());
@@ -108,9 +107,45 @@ ASSIMPReader::readVolumes(const std::string& source, double unit_length)  const
         vertices.emplace_back(Vertex(v[i].x*unit, v[i].y*unit, v[i].z*unit));
       }
       TessellatedSolid shape(name,vertices);
+      if ( name.empty() )  {
+        name = _toString(result.size(), "tessellated_%ld");
+      }
+
+      /// NOTE: IMPORTANT!
+      ///       ALWAYS add facets using the physical vertices!
+      ///       TGeoTessellated takes care that the vertex map is unique and
+      ///       assigns the proper indices to the facet.
       for(unsigned int i=0; i < mesh->mNumFaces; i++)  {
         const unsigned int* idx  = mesh->mFaces[i].mIndices;
-        shape->AddFacet(idx[0], idx[1], idx[2]);
+        bool degenerated = false;
+        if ( mesh->mFaces[i].mNumIndices == 3 )
+          degenerated = dd4hep::cad::facetIsDegenerated({vertices[idx[0]], vertices[idx[1]], vertices[idx[2]]});
+        else if ( mesh->mFaces[i].mNumIndices == 4 )
+          degenerated = dd4hep::cad::facetIsDegenerated({vertices[idx[0]], vertices[idx[1]], vertices[idx[2]], vertices[idx[3]]});
+
+        if ( degenerated )   {
+          printout(DEBUG, "ASSIMPReader", "+++ %s: Drop degenerated facet: %d %d %d",
+                   name.c_str(), idx[0], idx[1], idx[2]);
+        }
+#if ROOT_VERSION_CODE >= ROOT_VERSION(6,31,1)
+        else if ( mesh->mFaces[i].mNumIndices == 3 )   {
+          shape->AddFacet(vertices[idx[0]], vertices[idx[1]], vertices[idx[2]]);
+        }
+        else if ( mesh->mFaces[i].mNumIndices == 4 )   {
+          shape->AddFacet(vertices[idx[0]], vertices[idx[1]], vertices[idx[2]], vertices[idx[3]]);
+        }
+#else
+        else if ( mesh->mFaces[i].mNumIndices == 3 )   {
+          shape->AddFacet(idx[0], idx[1], idx[2]);
+        }
+        else if ( mesh->mFaces[i].mNumIndices == 4 )   {
+          shape->AddFacet(idx[0], idx[1], idx[2], idx[3]);
+        }
+#endif
+        else  {
+          printout(INFO, "ASSIMPReader", "+++ %s: Fancy facet with %d indices.",
+                   name.c_str(), mesh->mFaces[i].mNumIndices);
+        }
       }
       if ( shape->GetNfacets() > 2 )   {
         std::string mat_name;
@@ -124,14 +159,9 @@ ASSIMPReader::readVolumes(const std::string& source, double unit_length)  const
         if ( !mat.isValid() )   {
           printout(ERROR, "ASSIMPReader",
                    "+++ %s: No material named '%s' FOUND. Will use Air. [Missing material]",
-                   text, mat_name.c_str());
+                   name.c_str(), mat_name.c_str());
           mat = detector.air();
         }
-        if ( name.empty() )  {
-          ::snprintf(text,sizeof(text),"tessellated_%ld", result.size());
-          text[sizeof(text)-1] = 0;
-          name = text;
-        }
         Volume vol(name, Solid(shape.ptr()), mat);
         if ( mesh->HasVertexColors(0) )   {
           const aiColor4D* col = mesh->mColors[0];
@@ -146,6 +176,7 @@ ASSIMPReader::readVolumes(const std::string& source, double unit_length)  const
               }
             }
             if ( !vis.isValid() )   {
+              char text[1024];
               ::snprintf(text,sizeof(text),"vis_%s_%p", name.c_str(), (void*)vol.ptr());
               text[sizeof(text)-1] = 0;
               vis = VisAttr(text);
@@ -162,9 +193,9 @@ ASSIMPReader::readVolumes(const std::string& source, double unit_length)  const
         if ( dump_facets )   {
           for( size_t i=0, n=shape->GetNfacets(); i < n; ++i )   {
             const auto& facet = shape->GetFacet(i);
-            std::stringstream str = dd4hep::cad::streamFacet(facet, shape);
+            std::string str = dd4hep::cad::streamFacet(facet, shape);
             printout(ALWAYS,"ASSIMPReader","++ Facet %4ld : %s",
-                     i, str.str().c_str());
+                     i, str.c_str());
           }
         }
         result.emplace_back(std::unique_ptr<TGeoVolume>(vol.ptr()));

DDCAD/src/ASSIMPWriter.cpp

@@ -159,7 +159,7 @@ namespace  {
             continue;
           }
 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,31,1)
-          bool degenerated = dd4hep::cad::facetIsDegenerated(vertices);
+          bool degenerated = dd4hep::cad::facetIsDegenerated({vertices[vv0],vertices[vv1],vertices[vv2]});
 #else
           bool degenerated = true;
           TGeoFacet f(&vertices, 3, vv0, vv1, vv2);
@@ -426,8 +426,8 @@ int ASSIMPWriter::write(const std::string& file_name,
           ROOT::Geom::Vertex_t v1(vv[id[0]].x, vv[id[0]].y, vv[id[0]].z);
           ROOT::Geom::Vertex_t v2(vv[id[1]].x, vv[id[1]].y, vv[id[1]].z);
           ROOT::Geom::Vertex_t v3(vv[id[2]].x, vv[id[2]].y, vv[id[2]].z);
-          std::stringstream str = dd4hep::cad::streamVertices(v1, v2, v3);
-          printout(ALWAYS,"ASSIMPWriter","++ Facet %4ld : %s", j, str.str().c_str());
+          std::string str = dd4hep::cad::streamVertices(v1, v2, v3);
+          printout(ALWAYS,"ASSIMPWriter","++ Facet %4ld : %s", j, str.c_str());
         }
       }
       else   {

DDCore/include/DD4hep/Handle.h

@@ -239,10 +239,12 @@ namespace dd4hep {
 
   /// String conversions: boolean value to string  \ingroup DD4HEP_CORE
   std::string _toString(bool value);
-  /// String conversions: integer value to string  \ingroup DD4HEP_CORE
+  /// String conversions: short integer value to string  \ingroup DD4HEP_CORE
   std::string _toString(short value, const char* fmt = "%d");
   /// String conversions: integer value to string  \ingroup DD4HEP_CORE
   std::string _toString(int value, const char* fmt = "%d");
+  /// String conversions: unsigned long integer value to string  \ingroup DD4HEP_CORE
+  std::string _toString(unsigned long value, const char* fmt = "%ld");
   /// String conversions: float value to string  \ingroup DD4HEP_CORE
   std::string _toString(float value, const char* fmt = "%.17e");
   /// String conversions: double value to string  \ingroup DD4HEP_CORE

DDCore/src/Handle.cpp

@@ -344,6 +344,10 @@ namespace dd4hep  {
     return __to_string(value, fmt);
   }
 
+  string _toString(unsigned long value, const char* fmt) {
+    return __to_string(value, fmt);
+  }
+
   string _toString(float value, const char* fmt) {
     return __to_string(value, fmt);
   }

examples/ClientTests/compact/FiberTubeCalorimeter.xml

@@ -98,9 +98,9 @@
                 zmin="-world_side/2.+2*killthick+edgeoffset+DRcrystallength+EcalHcalgap"
                 z1="killthick"/>
       <structure>
-        <core1 name="core1" rmax="DRFiberFibwidthSCE" rmin="0.0" material="Polystyrene" vis="ScintVis" sensitive="yes"/> 
-        <core2 name="core2" rmax="DRFiberFibwidthSCE" rmin="0.0" material="Quartz" vis="CerenVis" sensitive="yes"/> 
-        <hole name="hole" rmax="(DRFiberFibwidthSCE+holeoverSCE)" rmin="0.0" material="Air" vis="holeVis" sensitive="yes"/> 
+        <core1  name="scintillator" rmax="DRFiberFibwidthSCE" rmin="0.0" material="Polystyrene" vis="ScintVis" sensitive="yes"/> 
+        <core2  name="quartz" rmax="DRFiberFibwidthSCE" rmin="0.0" material="Quartz" vis="CerenVis" sensitive="yes"/> 
+        <hole   name="hole" rmax="(DRFiberFibwidthSCE+holeoverSCE)" rmin="0.0" material="Air" vis="holeVis" sensitive="yes"/> 
         <absorb material="Brass" vis="AbsVis" sensitive="yes"/>
         <phdet1 name="phdet1" rmax="DRFiberFibwidthSCE" rmin="0.0" material="killMedia2" vis="phdetVis" sensitive="yes"/> 
         <phdet2 name="phdet2" rmax="DRFiberFibwidthSCE" rmin="0.0" material="killMedia3" vis="phdetVis" sensitive="yes"/>