From 349a6a436a75c5af9c64b4381c48dc04befc06e2 Mon Sep 17 00:00:00 2001
From: Markus Frank <Markus.Frank@cern.ch>
Date: Mon, 10 Apr 2017 18:37:35 +0200
Subject: [PATCH] Add isotope support when building elements

---
 DDCond/src/plugins/ConditionsUserPool.cpp     | 166 ++++++++++++++-
 DDCore/include/XML/UnicodeValues.h            |   1 +
 DDCore/src/LCDDImp.cpp                        |   4 +-
 DDCore/src/plugins/Compact2Objects.cpp        | 198 ++++++++++++++----
 .../SimpleDetector/compact/Simple_ILD.xml     |  12 ++
 examples/SimpleDetector/compact/elements.xml  |  87 ++++++--
 examples/SimpleDetector/compact/materials.xml |  66 ------
 7 files changed, 400 insertions(+), 134 deletions(-)

diff --git a/DDCond/src/plugins/ConditionsUserPool.cpp b/DDCond/src/plugins/ConditionsUserPool.cpp
index 69aa3de7a..08f4b1bd2 100644
--- a/DDCond/src/plugins/ConditionsUserPool.cpp
+++ b/DDCond/src/plugins/ConditionsUserPool.cpp
@@ -82,21 +82,26 @@ namespace DD4hep {
       virtual bool remove(const ConditionKey& key);
       /// Register a new condition to this pool
       virtual bool insert(Condition cond);
+
       /// Prepare user pool for usage (load, fill etc.) according to required IOV
       virtual Result prepare(const IOV&              required,
                              ConditionsSlice&        slice,
-                             void*                   user_param)  {
-        return prepare_VSN_1(required, slice, user_param);
-      }
-      /// Prepare user pool for usage (load, fill etc.) according to required IOV
-      Result prepare_VSN_1(const IOV&              required,
-                           ConditionsSlice&        slice,
-                           void*                   user_param);
+                             void*                   user_param);
+
       /// Evaluate and register all derived conditions from the dependency list
       virtual size_t compute(const Dependencies&     dependencies,
                              void*                   user_param,
                              bool                    force);
+
+      /// Prepare user pool for usage (load, fill etc.) according to required IOV
+      virtual Result load   (const IOV&              required,
+                             ConditionsSlice&        slice);
+      /// Prepare user pool for usage (load, fill etc.) according to required IOV
+      virtual Result compute(const IOV&              required,
+                             ConditionsSlice&        slice,
+                             void*                   user_param);
     };
+
   }    /* End namespace Conditions               */
 }      /* End namespace DD4hep                   */
 #endif /* DDCOND_CONDITIONSMAPPEDUSERPOOL_H      */
@@ -336,9 +341,9 @@ namespace {
 }
 
 template<typename MAPPING> UserPool::Result
-ConditionsMappedUserPool<MAPPING>::prepare_VSN_1(const IOV&              required, 
-                                                 ConditionsSlice&        slice,
-                                                 void*                   user_param)
+ConditionsMappedUserPool<MAPPING>::prepare(const IOV&              required, 
+                                           ConditionsSlice&        slice,
+                                           void*                   user_param)
 {
   typedef std::vector<pair<Condition::key_type,ConditionsDescriptor*> > _Missing;
   const auto& slice_cond = slice.conditions();
@@ -376,7 +381,7 @@ ConditionsMappedUserPool<MAPPING>::prepare_VSN_1(const IOV&              require
                                                 begin(m_conditions), end(m_conditions),
                                                 begin(calc_missing), COMP());
   int num_calc_miss = int(last_calc-begin(calc_missing));
-  printout(num_cond_miss==0 ? DEBUG : INFO,"UserPool",
+  printout(num_calc_miss==0 ? DEBUG : INFO,"UserPool",
            "Found %ld missing derived conditions out of %ld conditions.",
            num_calc_miss, m_conditions.size());
 
@@ -452,6 +457,145 @@ ConditionsMappedUserPool<MAPPING>::prepare_VSN_1(const IOV&              require
   return result;
 }
 
+template<typename MAPPING> UserPool::Result
+ConditionsMappedUserPool<MAPPING>::load(const IOV&              required, 
+                                        ConditionsSlice&        slice)
+{
+  typedef std::vector<pair<Condition::key_type,ConditionsDescriptor*> > _Missing;
+  const auto& slice_cond = slice.conditions();
+  auto&  slice_miss_cond = slice.missingConditions();
+  bool   do_load         = m_manager->doLoadConditions();
+  bool   do_output_miss  = m_manager->doOutputUnloaded();
+  IOV    pool_iov(required.iovType);
+  Result result;
+
+  // This is a critical operation, because we have to ensure the
+  // IOV pools are ONLY manipulated by the current thread.
+  // Otherwise the selection and the population are unsafe!
+  static mutex lock;
+  lock_guard<mutex> guard(lock);
+
+  m_conditions.clear();
+  slice_miss_cond.clear();
+  pool_iov.reset().invert();
+  m_iovPool->select(required, Operators::mapConditionsSelect(m_conditions), pool_iov);
+  m_iov = pool_iov;
+  _Missing cond_missing(slice_cond.size()+m_conditions.size());
+  _Missing::iterator last_cond = set_difference(begin(slice_cond),   end(slice_cond),
+                                                begin(m_conditions), end(m_conditions),
+                                                begin(cond_missing), COMP());
+  int num_cond_miss = int(last_cond-begin(cond_missing));
+  printout(num_cond_miss==0 ? DEBUG : INFO,"UserPool",
+           "Found %ld missing conditions out of %ld conditions.",
+           num_cond_miss, m_conditions.size());
+  result.loaded   = 0;
+  result.computed = 0;
+  result.missing  = num_cond_miss;
+  result.selected = m_conditions.size();
+  //
+  // Now we load the missing conditions from the conditions loader
+  //
+  if ( num_cond_miss > 0 )  {
+    if ( do_load )  {
+      ConditionsDataLoader::LoadedItems loaded;
+      size_t updates = m_loader->load_many(required, cond_missing, loaded, pool_iov);
+      if ( updates > 0 )  {
+        // Need to compute the intersection: All missing entries are required....
+        _Missing load_missing(cond_missing.size()+loaded.size());
+        // Note: cond_missing is already sorted (doc of 'set_difference'). No need to re-sort....
+        _Missing::iterator load_last = set_difference(begin(cond_missing), last_cond,
+                                                      begin(loaded), end(loaded),
+                                                      begin(load_missing), COMP());
+        int num_load_miss = int(load_last-begin(load_missing));
+        printout(num_load_miss==0 ? DEBUG : ERROR,"UserPool",
+                 "+++ %ld out of %d conditions CANNOT be loaded... [Not found by loader]",
+                 num_load_miss, int(loaded.size()));
+        if ( do_output_miss )  {
+          copy(begin(load_missing), load_last, inserter(slice_miss_cond, slice_miss_cond.begin()));
+        }
+        for_each(loaded.begin(),loaded.end(),Inserter<MAPPING>(m_conditions,&m_iov));
+        result.loaded  = slice_cond.size()-num_load_miss;
+        result.missing = num_load_miss;
+        if ( cond_missing.size() != loaded.size() )  {
+          // ERROR!
+        }
+      }
+    }
+    else if ( do_output_miss )  {
+      copy(begin(cond_missing), last_cond, inserter(slice_miss_cond, slice_miss_cond.begin()));
+    }
+  }
+  return result;
+}
+
+template<typename MAPPING> UserPool::Result
+ConditionsMappedUserPool<MAPPING>::compute(const IOV&              required, 
+                                           ConditionsSlice&        slice,
+                                           void*                   user_param)
+{
+  typedef std::vector<pair<Condition::key_type,ConditionsDescriptor*> > _Missing;
+  const auto& slice_calc = slice.derived();
+  auto&  slice_miss_calc = slice.missingDerivations();
+  bool   do_load         = m_manager->doLoadConditions();
+  bool   do_output_miss  = m_manager->doOutputUnloaded();
+  IOV    pool_iov(required.iovType);
+  Result result;
+
+  // This is a critical operation, because we have to ensure the
+  // IOV pools are ONLY manipulated by the current thread.
+  // Otherwise the selection and the population are unsafe!
+  static mutex lock;
+  lock_guard<mutex> guard(lock);
+
+  slice_miss_calc.clear();
+  _Missing calc_missing(slice_calc.size()+m_conditions.size());
+  _Missing::iterator last_calc = set_difference(begin(slice_calc),   end(slice_calc),
+                                                begin(m_conditions), end(m_conditions),
+                                                begin(calc_missing), COMP());
+  int num_calc_miss = int(last_calc-begin(calc_missing));
+  printout(num_calc_miss==0 ? DEBUG : INFO,"UserPool",
+           "Found %ld missing derived conditions out of %ld conditions.",
+           num_calc_miss, m_conditions.size());
+
+  result.loaded   = 0;
+  result.computed = 0;
+  result.missing  = num_calc_miss;
+  result.selected = m_conditions.size();
+  //
+  // Now we update the already existing dependencies, which have expired
+  //
+  if ( num_calc_miss > 0 )  {
+    if ( do_load )  {
+      ConditionsDependencyCollection deps(calc_missing.begin(), last_calc, _to_dep);
+      ConditionsDependencyHandler handler(m_manager, *this, deps, user_param);
+      for(auto i=begin(deps); i != end(deps); ++i)   {
+        const ConditionDependency* d = (*i).second.get();
+        typename MAPPING::iterator j = m_conditions.find(d->key());
+        // If we would know, that dependencies are only ONE level, we could skip this search....
+        if ( j == m_conditions.end() )  {
+          handler(d);
+          continue;
+        }
+        // printout(INFO,"UserPool","Already calcluated: %s",d->name());
+        continue;
+      }
+      result.computed = handler.num_callback;
+      result.missing -= handler.num_callback;
+      if ( do_output_miss && result.computed < deps.size() )  {
+        for(auto i=calc_missing.begin(); i != last_calc; ++i)   {
+          typename MAPPING::iterator j = m_conditions.find((*i).first);
+          if ( j == m_conditions.end() )
+            slice_miss_calc.insert(*i);
+        }
+      }
+    }
+    else if ( do_output_miss )  {
+      copy(begin(calc_missing), last_calc, inserter(slice_miss_calc, slice_miss_calc.begin()));
+    }
+  }
+  return result;
+}
+
 namespace {
   template <typename MAPPING>
   void* create_pool(Geometry::LCDD&, int argc, char** argv)  {
diff --git a/DDCore/include/XML/UnicodeValues.h b/DDCore/include/XML/UnicodeValues.h
index c1f009083..7df887717 100644
--- a/DDCore/include/XML/UnicodeValues.h
+++ b/DDCore/include/XML/UnicodeValues.h
@@ -86,6 +86,7 @@ UNICODE (cut);
 UNICODE (d);
 UNICODE (D);
 UNICODE (daughter);
+UNICODE (debug);
 UNICODE (define);
 UNICODE (delta);
 UNICODE (deltaphi);
diff --git a/DDCore/src/LCDDImp.cpp b/DDCore/src/LCDDImp.cpp
index ac31e6ef3..cf44afe28 100644
--- a/DDCore/src/LCDDImp.cpp
+++ b/DDCore/src/LCDDImp.cpp
@@ -115,14 +115,14 @@ LCDDImp::LCDDImp() : LCDDData(), LCDDLoad(this), m_buildType(BUILD_NONE)
   }
   {
     m_manager = gGeoManager;
-#if 0 //FIXME: eventually this should be set to 1 - needs fixes in examples ...
+    //#if 0 //FIXME: eventually this should be set to 1 - needs fixes in examples ...
     TGeoElementTable*	table = m_manager->GetElementTable();
     table->TGeoElementTable::~TGeoElementTable();
     new(table) TGeoElementTable();
     // This will initialize the table without filling:
     table->AddElement("VACUUM","VACUUM"   ,0,   0, 0.0);
     table->Print();
-#endif
+    //#endif
   }
   //if ( 0 == gGeoIdentity )
   {
diff --git a/DDCore/src/plugins/Compact2Objects.cpp b/DDCore/src/plugins/Compact2Objects.cpp
index 2a2e0823e..f0acb34f2 100644
--- a/DDCore/src/plugins/Compact2Objects.cpp
+++ b/DDCore/src/plugins/Compact2Objects.cpp
@@ -41,20 +41,24 @@ using namespace DD4hep::Geometry;
 
 namespace DD4hep {
   namespace Geometry {
-    struct Plugin;
-    struct Compact;
-    struct Includes;
-    struct GdmlFile;
-    struct Property;
-    struct XMLFile;
-    struct JsonFile;
-    struct AlignmentFile;
-    struct DetElementInclude {};
+    class Debug;
+    class Isotope;
+    class Plugin;
+    class Compact;
+    class Includes;
+    class GdmlFile;
+    class Property;
+    class XMLFile;
+    class JsonFile;
+    class AlignmentFile;
+    class DetElementInclude {};
   }
+  template <> void Converter<Debug>::operator()(xml_h element) const;
   template <> void Converter<Plugin>::operator()(xml_h element) const;
   template <> void Converter<Constant>::operator()(xml_h element) const;
   template <> void Converter<Material>::operator()(xml_h element) const;
   template <> void Converter<Atom>::operator()(xml_h element) const;
+  template <> void Converter<Isotope>::operator()(xml_h element) const;
   template <> void Converter<VisAttr>::operator()(xml_h element) const;
   template <> void Converter<AlignmentEntry>::operator()(xml_h element) const;
   template <> void Converter<Region>::operator()(xml_h element) const;
@@ -80,7 +84,14 @@ namespace {
     printout(ERROR, "Compact", msg.c_str());
     throw runtime_error(msg);
   }
-
+  bool s_debug_readout      = false;
+  bool s_debug_regions      = false;
+  bool s_debug_limits       = false;
+  bool s_debug_visattr      = false;
+  bool s_debug_isotopes     = false;
+  bool s_debug_elements     = false;
+  bool s_debug_materials    = false;
+  bool s_debug_segmentation = false;
 }
 
 static Ref_t create_ConstantField(lcdd_t& /* lcdd */, xml_h e) {
@@ -220,6 +231,23 @@ static long create_Compact(lcdd_t& lcdd, xml_h element) {
 }
 DECLARE_XML_DOC_READER(lccdd,create_Compact)
 
+/** Convert parser debug flags.
+ */
+template <> void Converter<Debug>::operator()(xml_h e) const {
+  for (xml_coll_t coll(e, _U(type)); coll; ++coll) {
+    string nam = coll.attr<string>(_U(name));
+    int    val = coll.attr<int>(_U(value));
+    if      ( nam.substr(0,6) == "isotop" ) s_debug_isotopes     = (0 != val);
+    else if ( nam.substr(0,6) == "elemen" ) s_debug_elements     = (0 != val);
+    else if ( nam.substr(0,6) == "materi" ) s_debug_materials    = (0 != val);
+    else if ( nam.substr(0,6) == "visatt" ) s_debug_visattr      = (0 != val);
+    else if ( nam.substr(0,6) == "region" ) s_debug_regions      = (0 != val);
+    else if ( nam.substr(0,6) == "readou" ) s_debug_readout      = (0 != val);
+    else if ( nam.substr(0,6) == "limits" ) s_debug_limits       = (0 != val);
+    else if ( nam.substr(0,6) == "segmen" ) s_debug_segmentation = (0 != val);
+  }
+}
+  
 /** Convert/execute plugin objects from the xml (plugins)
  *
  *
@@ -317,7 +345,8 @@ template <> void Converter<Material>::operator()(xml_h e) const {
       cout << " Density Value raw:" << dens_val << " normalized:" << (dens_val*dens_unit) << endl;
       dens_val *= dens_unit;
     }
-
+    printout(s_debug_materials ? ALWAYS : DEBUG, "Compact",
+             "++ Converting material %-16s  Density: %.3f.",matname, dens_val);
 #if 0
     cout << "Gev    " << XML::_toDouble(_Unicode(GeV)) << endl;
     cout << "sec    " << XML::_toDouble(_Unicode(second)) << endl;
@@ -411,31 +440,105 @@ template <> void Converter<Material>::operator()(xml_h e) const {
   }
 }
 
-/** Convert compact atom objects
+/** Convert compact isotope objects
  *
- *   <element Z="29" formula="Cu" name="Cu" >
+ *   <isotope name="C12" Z="2" N="12"/>
+ *     <atom unit="g/mole" value="12"/>
+ *   </isotope>
+ */
+template <> void Converter<Isotope>::operator()(xml_h e) const {
+  xml_dim_t isotope(e);
+  TGeoManager&      mgr  = lcdd.manager();
+  string            nam  = isotope.nameStr();
+  TGeoElementTable* tab  = mgr.GetElementTable();
+  TGeoIsotope*      iso  = tab->FindIsotope(nam.c_str());
+
+  // Create the isotope object in the event it is not yet present from the XML data
+  if ( !iso )   {
+    xml_ref_t atom(isotope.child(_U(atom)));
+    int    n     = isotope.attr<int>(_U(N));
+    int    z     = isotope.attr<int>(_U(Z));
+    double value = atom.attr<double>(_U(value));
+    string unit  = atom.attr<string>(_U(unit));
+    double a     = value * _multiply<double>(unit,"mol/g");
+    iso = new TGeoIsotope(nam.c_str(), z, n, a);
+    printout(s_debug_isotopes ? ALWAYS : DEBUG, "Compact",
+             "++ Converting isotope  %-16s  Z:%3d N:%3d A:%8.4f [g/mol]",
+             iso->GetName(), iso->GetZ(), iso->GetN(), iso->GetA());
+  }
+  else  {
+    printout(WARNING, "Compact",
+             "++ Isotope %-16s  Z:%3d N:%3d A:%8.4f [g/mol] ALREADY defined. [Ignore definition]",
+             iso->GetName(), iso->GetZ(), iso->GetN(), iso->GetA());
+  }
+}
+
+/** Convert compact atom objects (periodic elements)
+ *
+ *   <element Z="4" formula="Be" name="Be" >
+ *     <atom type="A" unit="g/mol" value="9.01218" />
+ *   </element>
+ *   or
+ *   <element name="C">
+ *     <fraction n="0.9893" ref="C12"/>
+ *     <fraction n="0.0107" ref="C13"/>
+ *   </element>
+ * 
+ *   Please note: 
+ *   Elements may consist of a mixture of isotopes!
  */
 template <> void Converter<Atom>::operator()(xml_h e) const {
-  xml_ref_t elem(e);
-  xml_tag_t eltname = elem.name();
-  TGeoManager& mgr = lcdd.manager();
-  TGeoElementTable* tab = mgr.GetElementTable();
-  TGeoElement* element = tab->FindElement(eltname.c_str());
-  if (!element) {
-    xml_ref_t atom(elem.child(_U(atom)));
-    string formula = elem.attr<string>(_U(formula));
-    double value   = atom.attr<double>(_U(value));
-    string unit    = atom.attr<string>(_U(unit));
-    int    z = elem.attr<int>(_U(Z));
-    double a = value*_multiply<double>(unit,"mol/g");
-    printout(DEBUG, "Compact", "++ Converting element  %-16s  [%-3s] Z:%3d A:%8.4f [g/mol]",
-             eltname.c_str(), formula.c_str(), z, a);
-    tab->AddElement(eltname.c_str(), formula.c_str(), z, a);
-    element = tab->FindElement(eltname.c_str());
-    if (!element) {
-      throw_print("Failed to properly insert the Element:" + eltname + " into the element table!");
+  xml_ref_t         elem(e);
+  xml_tag_t         name = elem.name();
+  TGeoManager&      mgr  = lcdd.manager();
+  TGeoElementTable* tab  = mgr.GetElementTable();
+  TGeoElement*      elt  = tab->FindElement(name.c_str());
+  if ( !elt ) {
+    if ( elem.hasChild(_U(atom)) )  {
+      xml_ref_t atom(elem.child(_U(atom)));
+      string formula = elem.attr<string>(_U(formula));
+      double value   = atom.attr<double>(_U(value));
+      string unit    = atom.attr<string>(_U(unit));
+      int    z       = elem.attr<int>(_U(Z));
+      double a       = value*_multiply<double>(unit,"mol/g");
+      printout(s_debug_elements ? ALWAYS : DEBUG, "Compact",
+               "++ Converting element  %-16s  [%-3s] Z:%3d A:%8.4f [g/mol]",
+               name.c_str(), formula.c_str(), z, a);
+      tab->AddElement(name.c_str(), formula.c_str(), z, a);
+    }
+    else  {
+      int num_isotopes = 0;
+      string formula   = elem.hasAttr(_U(formula)) ? elem.attr<string>(_U(formula)) : name.str();
+      for( xml_coll_t i(elem,_U(fraction)); i; ++i)
+        ++num_isotopes;
+      elt = new TGeoElement(name.c_str(), formula.c_str(), num_isotopes);
+      printout(s_debug_elements ? ALWAYS : DEBUG, "Compact",
+               "++ Converting element  %-16s  [%-3s] with %d isotopes.",
+               name.c_str(), formula.c_str(), num_isotopes);
+      tab->AddElement(elt);
+      for( xml_coll_t i(elem,_U(fraction)); i; ++i)  {
+        double frac = i.attr<double>(_U(n));
+        string ref  = i.attr<string>(_U(ref));
+        TGeoIsotope* iso = tab->FindIsotope(ref.c_str());
+        if ( !iso )  {
+          except("Compact","Element %s cannot be constructed. Isotope '%s' (fraction:%f) missing!",
+                 name.c_str(), ref.c_str(), frac);
+        }
+        printout(s_debug_elements ? ALWAYS : DEBUG, "Compact",
+                 "++ Converting element  %-16s  Add isotope: %-16s fraction:%.4f.",
+                 name.c_str(), ref.c_str(), frac);
+        elt->AddIsotope(iso, frac);
+      }
+    }
+    elt = tab->FindElement(name.c_str());
+    if (!elt) {
+      throw_print("Failed to properly insert the Element:"+name+" into the element table!");
     }
   }
+  else  {
+    printout(WARNING, "Compact", "++ Element %-16s  Z:%3d N:%3d A:%8.4f [g/mol] ALREADY defined. [Ignore definition]",
+             elt->GetName(), elt->Z(), elt->N(), elt->A());
+  }
 }
 
 /** Convert compact visualization attribute to LCDD visualization attribute
@@ -452,7 +555,9 @@ template <> void Converter<VisAttr>::operator()(xml_h e) const {
   float g = e.hasAttr(_U(g)) ? e.attr<float>(_U(g)) : 1.0f;
   float b = e.hasAttr(_U(b)) ? e.attr<float>(_U(b)) : 1.0f;
 
-  printout(DEBUG, "Compact", "++ Converting VisAttr  structure: %s.",attr.name());
+  printout(s_debug_visattr ? ALWAYS : DEBUG, "Compact",
+           "++ Converting VisAttr  structure: %-16s. R=%.3f G=%.3f B=%.3f",
+           attr.name(), r, g, b);
   attr.setColor(r, g, b);
   if (e.hasAttr(_U(alpha)))
     attr.setAlpha(e.attr<float>(_U(alpha)));
@@ -528,7 +633,8 @@ template <> void Converter<Region>::operator()(xml_h elt) const {
   xml_attr_t store_secondaries = elt.attr_nothrow(_U(store_secondaries));
   double ene = e.eunit(1.0), len = e.lunit(1.0);
 
-  printout(DEBUG, "Compact", "++ Converting region   structure: %s.",region.name());
+  printout(s_debug_regions ? ALWAYS : DEBUG, "Compact",
+           "++ Converting region   structure: %s.",region.name());
   if ( cut )  {
     region.setCut(elt.attr<double>(cut)*len);
   }
@@ -561,6 +667,8 @@ template <> void Converter<Segmentation>::operator()(xml_h seg) const {
   if (segment.isValid()) {
     typedef Segmentation::Parameters _PARS;
     const _PARS& pars = segment.parameters();
+    printout(s_debug_segmentation ? ALWAYS : DEBUG, "Compact",
+             "++ Converting segmentation structure: %s of type %s.",name.c_str(),type.c_str());
     for(_PARS::const_iterator it = pars.begin(); it != pars.end(); ++it) {
       Segmentation::Parameter p = *it;
       XML::Strng_t pNam(p->name());
@@ -575,7 +683,8 @@ template <> void Converter<Segmentation>::operator()(xml_h seg) const {
         } else if ( pType.compare("doublevec") == 0 ) {
           vector<double> valueVector;
           string par = seg.attr<string>(pNam);
-          printout(DEBUG, "Compact", "++ Converting this string structure: %s.",par.c_str());
+          printout(s_debug_segmentation ? ALWAYS : DEBUG, "Compact",
+                   "++ Converting this string structure: %s.",par.c_str());
           vector<string> elts = DDSegmentation::splitString(par);
           for (vector<string>::const_iterator j = elts.begin(); j != elts.end(); ++j) {
             if ((*j).empty()) continue;
@@ -615,7 +724,8 @@ template <> void Converter<Segmentation>::operator()(xml_h seg) const {
           XML::dump_tree(sub,tree);
           throw_print("Nested segmentations: Invalid key specification:"+tree.str());
         }
-        printout(DEBUG,"Compact","++ Segmentation [%s/%s]: Add sub-segmentation %s [%s]",
+        printout(s_debug_segmentation ? ALWAYS : DEBUG,"Compact",
+                 "++ Segmentation [%s/%s]: Add sub-segmentation %s [%s]",
                  name.c_str(), type.c_str(), 
                  sub_seg->segmentation->name().c_str(),
                  sub_seg->segmentation->type().c_str());
@@ -642,8 +752,6 @@ template <> void Converter<Readout>::operator()(xml_h e) const {
   std::pair<Segmentation,IDDescriptor> opt;
   Readout ro(name);
   
-  printout(DEBUG, "Compact", "++ Converting readout  structure: %s.",ro.name());
-  
   if (id) {
     //  <id>system:6,barrel:3,module:4,layer:8,slice:5,x:32:-16,y:-16</id>
     opt.second = IDDescriptor(id.text());
@@ -664,6 +772,10 @@ template <> void Converter<Readout>::operator()(xml_h e) const {
     ro.setIDDescriptor(opt.second);
   }
   
+  printout(s_debug_readout ? ALWAYS : DEBUG,
+           "Compact", "++ Converting readout  structure: %-16s. %s%s",
+           ro.name(), id ? "ID: " : "", id ? id.text().c_str() : "");
+  
   for(xml_coll_t colls(e,_U(hits_collections)); colls; ++colls)   {
     string hits_key;
     if ( colls.hasAttr(_U(key)) ) hits_key = colls.attr<string>(_U(key));
@@ -688,7 +800,8 @@ template <> void Converter<Readout>::operator()(xml_h e) const {
         XML::dump_tree(e,tree);
         throw_print("Reaout: Invalid specificatrion for multiple hit collections."+tree.str());
       }
-      printout(DEBUG,"Compact","++ Readout[%s]: Add hit collection %s [%s]  %d-%d",
+      printout(s_debug_readout ? ALWAYS : DEBUG,"Compact",
+               "++ Readout[%s]: Add hit collection %s [%s]  %d-%d",
                ro.name(), coll_name.c_str(), coll_key.c_str(), key_min, key_max);
       HitCollection hits(coll_name, coll_key, key_min, key_max);
       ro->hits.push_back(hits);
@@ -705,7 +818,8 @@ template <> void Converter<Readout>::operator()(xml_h e) const {
  */
 template <> void Converter<LimitSet>::operator()(xml_h e) const {
   LimitSet ls(e.attr<string>(_U(name)));
-  printout(DEBUG, "Compact", "++ Converting LimitSet structure: %s.",ls.name());
+  printout(s_debug_limits ? ALWAYS : DEBUG, "Compact",
+           "++ Converting LimitSet structure: %s.",ls.name());
   for (xml_coll_t c(e, _U(limit)); c; ++c) {
     Limit limit;
     limit.particles = c.attr<string>(_U(particles));
@@ -953,9 +1067,10 @@ template <> void Converter<DetElement>::operator()(xml_h element) const {
 }
 
 /// Read material entries from a seperate file in one of the include sections of the geometry
-template <> void Converter<GdmlFile>::operator()(xml_h element) const {
+template <> void Converter<GdmlFile>::operator()(xml_h element) const   {
   XML::DocumentHolder doc(XML::DocumentHandler().load(element, element.attr_value(_U(ref))));
   xml_h materials = doc.root();
+  xml_coll_t(materials, _U(isotope)).for_each(Converter<Isotope>(this->lcdd,0,0));
   xml_coll_t(materials, _U(element)).for_each(Converter<Atom>(this->lcdd));
   xml_coll_t(materials, _U(material)).for_each(Converter<Material>(this->lcdd));
 }
@@ -1029,6 +1144,9 @@ template <> void Converter<Compact>::operator()(xml_h element) const {
   bool open_geometry  = true;
   bool close_geometry = true;
 
+  if (element.hasChild(_U(debug)))
+    (Converter<Debug>(lcdd))(xml_h(compact.child(_U(debug))));
+  
   if ( steer_geometry )   {
     xml_elt_t steer = compact.child(_U(geometry));
     if ( steer.hasAttr(_U(open))  ) open_geometry  = steer.attr<bool>(_U(open));
diff --git a/examples/SimpleDetector/compact/Simple_ILD.xml b/examples/SimpleDetector/compact/Simple_ILD.xml
index 32ce0dc6e..eb7913d49 100644
--- a/examples/SimpleDetector/compact/Simple_ILD.xml
+++ b/examples/SimpleDetector/compact/Simple_ILD.xml
@@ -13,6 +13,18 @@
            so far only VXD and (pixel) SIT ...
 	</comment>        
     </info>
+
+    <!-- Debug flags for processing the compact xml file  -->
+    <debug>
+      <type name="isotopes"     value="0"/>
+      <type name="elements"     value="0"/>
+      <type name="materials"    value="0"/>
+      <type name="visattr"      value="0"/>
+      <type name="regions"      value="0"/>
+      <type name="readout"      value="0"/>
+      <type name="limits"       value="0"/>
+      <type name="segmentation" value="0"/>
+    </debug>
  
     <includes>
         <gdmlFile  ref="elements.xml"/>
diff --git a/examples/SimpleDetector/compact/elements.xml b/examples/SimpleDetector/compact/elements.xml
index e714c3a5c..9baedbfd5 100644
--- a/examples/SimpleDetector/compact/elements.xml
+++ b/examples/SimpleDetector/compact/elements.xml
@@ -1,4 +1,76 @@
 <materials>
+
+    <isotope N="1" Z="1" name="H1">
+      <atom unit="g/mole" value="1.00782503081372"/>
+    </isotope>
+    <isotope N="2" Z="1" name="H2">
+      <atom unit="g/mole" value="2.01410199966617"/>
+    </isotope>
+    <element name="H">
+      <fraction n="0.999885" ref="H1"/>
+      <fraction n="0.000115" ref="H2"/>
+    </element>
+
+    <isotope N="12" Z="6" name="C12">
+      <atom unit="g/mole" value="12"/>
+    </isotope>
+    <isotope N="13" Z="6" name="C13">
+      <atom unit="g/mole" value="13.0034"/>
+    </isotope>
+    <element name="C">
+      <fraction n="0.9893" ref="C12"/>
+      <fraction n="0.0107" ref="C13"/>
+    </element>
+
+    <isotope N="14" Z="7" name="N14">
+      <atom unit="g/mole" value="14.0031"/>
+    </isotope>
+    <isotope N="15" Z="7" name="N15">
+      <atom unit="g/mole" value="15.0001"/>
+    </isotope>
+    <element name="N">
+      <fraction n="0.99632" ref="N14"/>
+      <fraction n="0.00368" ref="N15"/>
+    </element>
+
+    <isotope N="16" Z="8" name="O16">
+      <atom unit="g/mole" value="15.9949"/>
+    </isotope>
+    <isotope N="17" Z="8" name="O17">
+      <atom unit="g/mole" value="16.9991"/>
+    </isotope>
+    <isotope N="18" Z="8" name="O18">
+      <atom unit="g/mole" value="17.9992"/>
+    </isotope>
+    <element name="O">
+      <fraction n="0.99757" ref="O16"/>
+      <fraction n="0.00038" ref="O17"/>
+      <fraction n="0.00205" ref="O18"/>
+    </element>
+
+    <isotope N="46" Z="22" name="Ti46">
+      <atom unit="g/mole" value="45.9526"/>
+    </isotope>
+    <isotope N="47" Z="22" name="Ti47">
+      <atom unit="g/mole" value="46.9518"/>
+    </isotope>
+    <isotope N="48" Z="22" name="Ti48">
+      <atom unit="g/mole" value="47.9479"/>
+    </isotope>
+    <isotope N="49" Z="22" name="Ti49">
+      <atom unit="g/mole" value="48.9479"/>
+    </isotope>
+    <isotope N="50" Z="22" name="Ti50">
+      <atom unit="g/mole" value="49.9448"/>
+    </isotope>
+    <element name="Ti">
+      <fraction n="0.0825" ref="Ti46"/>
+      <fraction n="0.0744" ref="Ti47"/>
+      <fraction n="0.7372" ref="Ti48"/>
+      <fraction n="0.0541" ref="Ti49"/>
+      <fraction n="0.0518" ref="Ti50"/>
+    </element>
+
  <element Z="89" formula="Ac" name="Ac" >
   <atom type="A" unit="g/mol" value="227.028" />
  </element>
@@ -125,9 +197,6 @@
   <D type="density" unit="g/cm3" value="0.0070721" />
   <composite n="1" ref="Br" />
  </material>
- <element Z="6" formula="C" name="C" >
-  <atom type="A" unit="g/mol" value="12.0107" />
- </element>
  <material formula="C" name="Carbon" state="solid" >
   <RL type="X0" unit="cm" value="21.3485" />
   <NIL type="lambda" unit="cm" value="40.1008" />
@@ -305,9 +374,6 @@
   <D type="density" unit="g/cm3" value="5.323" />
   <composite n="1" ref="Ge" />
  </material>
- <element Z="1" formula="H" name="H" >
-  <atom type="A" unit="g/mol" value="1.00794" />
- </element>
  <material formula="H" name="Hydrogen" state="gas" >
   <RL type="X0" unit="cm" value="752776" />
   <NIL type="lambda" unit="cm" value="421239" />
@@ -449,9 +515,6 @@
   <D type="density" unit="g/cm3" value="10.22" />
   <composite n="1" ref="Mo" />
  </material>
- <element Z="7" formula="N" name="N" >
-  <atom type="A" unit="g/mol" value="14.0068" />
- </element>
  <material formula="N" name="Nitrogen" state="gas" >
   <RL type="X0" unit="cm" value="32602.2" />
   <NIL type="lambda" unit="cm" value="72430.3" />
@@ -512,9 +575,6 @@
   <D type="density" unit="g/cm3" value="20.25" />
   <composite n="1" ref="Np" />
  </material>
- <element Z="8" formula="O" name="O" >
-  <atom type="A" unit="g/mol" value="15.9994" />
- </element>
  <material formula="O" name="Oxygen" state="gas" >
   <RL type="X0" unit="cm" value="25713.8" />
   <NIL type="lambda" unit="cm" value="66233.9" />
@@ -782,9 +842,6 @@
   <D type="density" unit="g/cm3" value="11.72" />
   <composite n="1" ref="Th" />
  </material>
- <element Z="22" formula="Ti" name="Ti" >
-  <atom type="A" unit="g/mol" value="47.8667" />
- </element>
  <material formula="Ti" name="Titanium" state="solid" >
   <RL type="X0" unit="cm" value="3.5602" />
   <NIL type="lambda" unit="cm" value="27.9395" />
diff --git a/examples/SimpleDetector/compact/materials.xml b/examples/SimpleDetector/compact/materials.xml
index c5fd302b5..200798e1a 100644
--- a/examples/SimpleDetector/compact/materials.xml
+++ b/examples/SimpleDetector/compact/materials.xml
@@ -1,48 +1,4 @@
   <materials>
-    <isotope N="1" Z="1" name="H1">
-      <atom unit="g/mole" value="1.00782503081372"/>
-    </isotope>
-    <isotope N="2" Z="1" name="H2">
-      <atom unit="g/mole" value="2.01410199966617"/>
-    </isotope>
-    <element name="H">
-      <fraction n="0.999885" ref="H1"/>
-      <fraction n="0.000115" ref="H2"/>
-    </element>
-    <isotope N="12" Z="6" name="C12">
-      <atom unit="g/mole" value="12"/>
-    </isotope>
-    <isotope N="13" Z="6" name="C13">
-      <atom unit="g/mole" value="13.0034"/>
-    </isotope>
-    <element name="C">
-      <fraction n="0.9893" ref="C12"/>
-      <fraction n="0.0107" ref="C13"/>
-    </element>
-    <isotope N="14" Z="7" name="N14">
-      <atom unit="g/mole" value="14.0031"/>
-    </isotope>
-    <isotope N="15" Z="7" name="N15">
-      <atom unit="g/mole" value="15.0001"/>
-    </isotope>
-    <element name="N">
-      <fraction n="0.99632" ref="N14"/>
-      <fraction n="0.00368" ref="N15"/>
-    </element>
-    <isotope N="16" Z="8" name="O16">
-      <atom unit="g/mole" value="15.9949"/>
-    </isotope>
-    <isotope N="17" Z="8" name="O17">
-      <atom unit="g/mole" value="16.9991"/>
-    </isotope>
-    <isotope N="18" Z="8" name="O18">
-      <atom unit="g/mole" value="17.9992"/>
-    </isotope>
-    <element name="O">
-      <fraction n="0.99757" ref="O16"/>
-      <fraction n="0.00038" ref="O17"/>
-      <fraction n="0.00205" ref="O18"/>
-    </element>
     <material name="G4_KAPTON" state="solid">
       <MEE unit="eV" value="79.6"/>
       <D unit="g/cm3" value="1.42"/>
@@ -103,28 +59,6 @@
       <atom unit="g/mole" value="9.01218"/>
       <fraction n="1" ref="Be"/>
     </material>
-    <isotope N="46" Z="22" name="Ti46">
-      <atom unit="g/mole" value="45.9526"/>
-    </isotope>
-    <isotope N="47" Z="22" name="Ti47">
-      <atom unit="g/mole" value="46.9518"/>
-    </isotope>
-    <isotope N="48" Z="22" name="Ti48">
-      <atom unit="g/mole" value="47.9479"/>
-    </isotope>
-    <isotope N="49" Z="22" name="Ti49">
-      <atom unit="g/mole" value="48.9479"/>
-    </isotope>
-    <isotope N="50" Z="22" name="Ti50">
-      <atom unit="g/mole" value="49.9448"/>
-    </isotope>
-    <element name="Ti">
-      <fraction n="0.0825" ref="Ti46"/>
-      <fraction n="0.0744" ref="Ti47"/>
-      <fraction n="0.7372" ref="Ti48"/>
-      <fraction n="0.0541" ref="Ti49"/>
-      <fraction n="0.0518" ref="Ti50"/>
-    </element>
     <material name="G4_Ti" state="solid">
       <MEE unit="eV" value="233"/>
       <D unit="g/cm3" value="4.54"/>
-- 
GitLab