Geant4Converter.cpp

        Geant4AssemblyVolume::Chain chain;
        chain.emplace_back(node);
        ass->imprint(*this, node, chain, ass, (*volIt).second, transform, copy, checkOverlaps);
        return nullptr;
      }
      else if ( node != info.manager->GetTopNode() && volIt == info.g4Volumes.end() )  {
        //throw std::logic_error("Geant4Converter: Invalid mother volume found!");
      }
      PlacedVolume pv(node);
      const auto*  pv_data = pv.data();
      G4LogicalVolume* g4vol = info.g4Volumes[vol];
      G4LogicalVolume* g4mot = info.g4Volumes[mot_vol];
      G4PhysicalVolumesPair pvPlaced  { nullptr, nullptr };

      if ( pv_data && pv_data->params && (pv_data->params->flags&Volume::REPLICATED) )   {
        EAxis  axis = kUndefined;
        double width = 0e0, offset = 0e0;
        auto flags = pv_data->params->flags;
        auto count = pv_data->params->trafo1D.second;
        auto start = pv_data->params->start.Translation().Vect();
        auto delta = pv_data->params->trafo1D.first.Translation().Vect();

        if ( flags&Volume::X_axis )
        { axis = kXAxis; width = delta.X(); offset = start.X(); }
        else if ( flags&Volume::Y_axis )
        { axis = kYAxis; width = delta.Y(); offset = start.Y(); }
        else if ( flags&Volume::Z_axis )
        { axis = kZAxis; width = delta.Z(); offset = start.Z(); }
        else
          except("Geant4Converter",
                 "++ Replication around unknown axis is not implemented. flags: %16X", flags);
        printout(INFO,"Geant4Converter","++ Replicate: Axis: %ld Count: %ld offset: %f width: %f",
                 axis, count, offset, width);
        auto* g4pv = new G4PVReplica(name,      // its name
                                     g4vol,     // its logical volume
                                     g4mot,     // its mother (logical) volume
                                     axis,      // its replication axis
                                     count,     // Number of replicas
                                     width,     // Distance between 2 replicas
                                     offset);   // Placement offset in axis direction
        pvPlaced = { g4pv, nullptr };
#if 0
        pvPlaced =
          G4ReflectionFactory::Instance()->Replicate(name,      // its name
                                                     g4vol,     // its logical volume
                                                     g4mot,     // its mother (logical) volume
                                                     axis,      // its replication axis
                                                     count,     // Number of replicas
                                                     width,     // Distance between 2 replicas
                                                     offset);   // Placement offset in axis direction
        /// Update replica list to avoid additional conversions...
        auto* g4pv = pvPlaced.second ? pvPlaced.second : pvPlaced.first;
#endif
        for( auto& handle : pv_data->params->placements )
          info.g4Placements[handle.ptr()] = g4pv;
      }
      else if ( pv_data && pv_data->params )   {
        auto*  g4par = new Geant4PlacementParameterisation(pv);
        auto*  g4pv  = new G4PVParameterised(name,              // its name
                                             g4vol,             // its logical volume
                                             g4mot,             // its mother (logical) volume
                                             g4par->axis(),     // its replication axis
                                             g4par->count(),    // Number of replicas
                                             g4par);            // G4 parametrization
        pvPlaced = { g4pv, nullptr };
        /// Update replica list to avoid additional conversions...
        for( auto& handle : pv_data->params->placements )
          info.g4Placements[handle.ptr()] = g4pv;
      }
      else    {
        pvPlaced =
          G4ReflectionFactory::Instance()->Place(transform,     // no rotation
                                                 name,          // its name
                                                 g4vol,         // its logical volume
                                                 g4mot,         // its mother (logical) volume
                                                 false,         // no boolean operations
                                                 copy,          // its copy number
                                                 checkOverlaps);
      }
      printout(debugReflections||debugPlacements ? ALWAYS : lvl, "Geant4Converter",
               "++ Place %svolume %-12s in mother %-12s "
               "Tr:x=%8.1f y=%8.1f z=%8.1f   Scale:x=%4.2f y=%4.2f z=%4.2f",
               node_is_reflected ? "REFLECTED " : "", _v.name(),
               mot_vol ? mot_vol->GetName() : "<unknown>",
               transform.dx(), transform.dy(), transform.dz(),
               scale.xx(), scale.yy(), scale.zz());
      // First 2 cases can be combined.
      // Leave them separated for debugging G4ReflectionFactory for now...
      if ( node_is_reflected  && !pvPlaced.second )
        return info.g4Placements[node] = pvPlaced.first;
      else if ( !node_is_reflected && !pvPlaced.second )
        return info.g4Placements[node] = pvPlaced.first;
      // Now deal with valid pvPlaced.second ...
      if ( node_is_reflected )
        return info.g4Placements[node] = pvPlaced.first;
      else if ( !node_is_reflected )
        return info.g4Placements[node] = pvPlaced.first;
      g4 = pvPlaced.second ? pvPlaced.second : pvPlaced.first;
    }
    info.g4Placements[node] = g4;
    printout(ERROR, "Geant4Converter", "++ DEAD code. Should not end up here!");
  }
  return g4;
}

/// Convert the geometry type region into the corresponding Geant4 object(s).
void* Geant4Converter::handleRegion(Region region, const std::set<const TGeoVolume*>& /* volumes */) const {
  G4Region* g4 = data().g4Regions[region];
  if ( !g4 ) {
    PrintLevel lvl = debugRegions ? ALWAYS : outputLevel;
    Region r = region;
    g4 = new G4Region(region.name());

    // create region info with storeSecondaries flag
    if( not r.wasThresholdSet() and r.storeSecondaries() ) {
      throw std::runtime_error("G4Region: StoreSecondaries is True, but no explicit threshold set:");
    }
    printout(lvl, "Geant4Converter", "++ Setting up region: %s", r.name());
    G4UserRegionInformation* info = new G4UserRegionInformation();
    info->region = r;
    info->threshold = r.threshold()*CLHEP::MeV/units::MeV;
    info->storeSecondaries = r.storeSecondaries();
    g4->SetUserInformation(info);

    printout(lvl, "Geant4Converter", "++ Converted region settings of:%s.", r.name());
    std::vector < std::string > &limits = r.limits();
    G4ProductionCuts* cuts = 0;
    // set production cut
    if( not r.useDefaultCut() ) {
      cuts = new G4ProductionCuts();
      cuts->SetProductionCut(r.cut()*CLHEP::mm/units::mm);
      printout(lvl, "Geant4Converter", "++ %s: Using default cut: %f [mm]",
               r.name(), r.cut()*CLHEP::mm/units::mm);
    }
    for( const auto& nam : limits )  {
      LimitSet ls = m_detDesc.limitSet(nam);
      if ( ls.isValid() ) {
        const LimitSet::Set& cts = ls.cuts();
        for (const auto& c : cts )   {
          int pid = 0;
          if ( c.particles == "*" ) pid = -1;
          else if ( c.particles == "e-"     ) pid = idxG4ElectronCut;
          else if ( c.particles == "e+"     ) pid = idxG4PositronCut;
          else if ( c.particles == "e[+-]"  ) pid = -idxG4PositronCut-idxG4ElectronCut;
          else if ( c.particles == "e[-+]"  ) pid = -idxG4PositronCut-idxG4ElectronCut;
          else if ( c.particles == "gamma"  ) pid = idxG4GammaCut;
          else if ( c.particles == "proton" ) pid = idxG4ProtonCut;
          else throw std::runtime_error("G4Region: Invalid production cut particle-type:" + c.particles);
          if ( !cuts ) cuts = new G4ProductionCuts();
          if ( pid == -(idxG4PositronCut+idxG4ElectronCut) )  {
            cuts->SetProductionCut(c.value*CLHEP::mm/units::mm, idxG4PositronCut);
            cuts->SetProductionCut(c.value*CLHEP::mm/units::mm, idxG4ElectronCut);
          }
          else  {
            cuts->SetProductionCut(c.value*CLHEP::mm/units::mm, pid);
          }
          printout(lvl, "Geant4Converter", "++ %s: Set cut  [%s/%d] = %f [mm]",
                   r.name(), c.particles.c_str(), pid, c.value*CLHEP::mm/units::mm);
        }
        bool found = false;
        const auto& lm = data().g4Limits;
        for (const auto& j : lm )   {
          if (nam == j.first->GetName()) {
            g4->SetUserLimits(j.second);
            printout(lvl, "Geant4Converter", "++ %s: Set limits %s to region type %s",
                     r.name(), nam.c_str(), j.second->GetType().c_str());
            found = true;
            break;
          }
        }
        if ( found )   {
          continue;
        }
      }
      except("Geant4Converter", "++ G4Region: Failed to resolve limitset: " + nam);
    }
    /// Assign cuts to region if they were created
    if ( cuts ) g4->SetProductionCuts(cuts);
    data().g4Regions[region] = g4;
  }
  return g4;
}

/// Convert the geometry type LimitSet into the corresponding Geant4 object(s).
void* Geant4Converter::handleLimitSet(LimitSet limitset, const std::set<const TGeoVolume*>& /* volumes */) const {
  G4UserLimits* g4 = data().g4Limits[limitset];
  if ( !g4 ) {
    PrintLevel lvl = debugLimits || debugRegions ? ALWAYS : outputLevel;
    struct LimitPrint  {
      const LimitSet& ls;
      LimitPrint(const LimitSet& lset) : ls(lset) {}
      const LimitPrint& operator()(const std::string& pref, const Geant4UserLimits::Handler& h)  const {
        if ( !h.particleLimits.empty() )  {
          printout(ALWAYS,"Geant4Converter",
                   "+++ LimitSet: Explicit Limit %s.%s applied for particles:",ls.name(), pref.c_str());
          for(const auto& p : h.particleLimits)
            printout(ALWAYS,"Geant4Converter","+++ LimitSet:    Particle type: %-18s PDG: %-6d : %f",
                     p.first->GetParticleName().c_str(), p.first->GetPDGEncoding(), p.second);
        }
        else if ( h.defaultValue > std::numeric_limits<double>::epsilon() )  {
          printout(ALWAYS,"Geant4Converter",
                   "+++ LimitSet: Implicit Limit %s.%s for wildcard particles: %f",
                   ls.name(), pref.c_str(), float(h.defaultValue));
        }
        return *this;
      }
    };
    Geant4UserLimits* limits = new Geant4UserLimits(limitset);
    g4 = limits;
    printout(lvl, "Geant4Converter",
             "++ Successfully converted LimitSet: %s [%ld cuts, %ld limits]",
             limitset.name(), limitset.cuts().size(), limitset.limits().size());
    if ( debugRegions || debugLimits )    {
      LimitPrint print(limitset);
      print("maxTime",    limits->maxTime)
        ("minEKine",      limits->minEKine)
        ("minRange",      limits->minRange)
        ("maxStepLength", limits->maxStepLength)
        ("maxTrackLength",limits->maxTrackLength);
    }
    data().g4Limits[limitset] = g4;
  }
  return g4;
}

/// Convert the geometry visualisation attributes to the corresponding Geant4 object(s).
void* Geant4Converter::handleVis(const std::string& /* name */, VisAttr attr) const {
  Geant4GeometryInfo& info = data();
  G4VisAttributes*    g4   = info.g4Vis[attr];
  if ( !g4 ) {
    float red = 0, green = 0, blue = 0;
    int   style = attr.lineStyle();
    attr.rgb(red, green, blue);
    g4 = new G4VisAttributes(attr.visible(), G4Colour(red, green, blue, attr.alpha()));
    //g4->SetLineWidth(attr->GetLineWidth());
    g4->SetDaughtersInvisible(!attr.showDaughters());
    if ( style == VisAttr::SOLID ) {
      g4->SetLineStyle(G4VisAttributes::unbroken);
      g4->SetForceWireframe(false);
      g4->SetForceSolid(true);
    }
    else if ( style == VisAttr::WIREFRAME || style == VisAttr::DASHED ) {
      g4->SetLineStyle(G4VisAttributes::dashed);
      g4->SetForceSolid(false);
      g4->SetForceWireframe(true);
    }
    info.g4Vis[attr] = g4;
  }
  return g4;
}

/// Handle the geant 4 specific properties
void Geant4Converter::handleProperties(Detector::Properties& prp) const {
  std::map < std::string, std::string > processors;
  static int s_idd = 9999999;
  for( const auto& [nam, vals] : prp ) {
    if ( nam.substr(0, 6) == "geant4" ) {
      auto id_it = vals.find("id");
      std::string id = (id_it == vals.end()) ? _toString(++s_idd,"%d") : (*id_it).second;
      processors.emplace(id, nam);
    }
  }
  for( const auto& p : processors ) {
    const GeoHandler* hdlr = this;
    const Detector::PropertyValues& vals = prp[p.second];
    auto iter = vals.find("type");
    if ( iter != vals.end() )  {
      std::string type = iter->second;
      std::string tag  = type + "_Geant4_action";
      Detector* det = const_cast<Detector*>(&m_detDesc);
      long      res = PluginService::Create<long>(tag, det, hdlr, &vals);
      if ( 0 == res ) {
        throw std::runtime_error("Failed to locate plugin to interprete files of type"
                                 " \"" + tag + "\" - no factory:" + type);
      }
      res = *(long*)res;
      if ( res != 1 ) {
        throw std::runtime_error("Failed to invoke the plugin " + tag + " of type " + type);
      }
      printout(outputLevel, "Geant4Converter", "+++++ Executed Successfully Geant4 setup module *%s*.", type.c_str());
      continue;
    }
    printout(outputLevel, "Geant4Converter", "+++++ FAILED to execute Geant4 setup module *%s*.", p.second.c_str());    
  }
}

/// Convert the geometry type material into the corresponding Geant4 object(s).
void* Geant4Converter::handleMaterialProperties(TObject* mtx) const    {
  Geant4GeometryInfo& info   = data();
  TGDMLMatrix*        matrix = (TGDMLMatrix*)mtx;
  const char*         cptr   = ::strstr(matrix->GetName(), GEANT4_TAG_IGNORE);
  Geant4GeometryInfo::PropertyVector* g4 = info.g4OpticalProperties[matrix];

  if ( nullptr != cptr )   {  // Check if the property should not be passed to Geant4
    printout(INFO,"Geant4MaterialProperties","++ Ignore property %s [%s].",
             matrix->GetName(), matrix->GetTitle());	     
    return nullptr;
  }
  cptr = ::strstr(matrix->GetTitle(), GEANT4_TAG_IGNORE);
  if ( nullptr != cptr )   {  // Check if the property should not be passed to Geant4
    printout(INFO,"Geant4MaterialProperties","++ Ignore property %s [%s].",
             matrix->GetName(), matrix->GetTitle());
    return nullptr;
  }
  
  if ( !g4 )  {
    PrintLevel lvl = debugMaterials ? ALWAYS : outputLevel;
    g4 = new Geant4GeometryInfo::PropertyVector();
    std::size_t rows = matrix->GetRows();
    g4->name    = matrix->GetName();
    g4->title   = matrix->GetTitle();
    g4->bins.reserve(rows);
    g4->values.reserve(rows);
    for( std::size_t i=0; i<rows; ++i )   {
      g4->bins.emplace_back(matrix->Get(i,0)  /*   *CLHEP::eV/units::eV   */);
      g4->values.emplace_back(matrix->Get(i,1));
    }
    printout(lvl, "Geant4Converter",
             "++ Successfully converted material property:%s : %s [%ld rows]",
             matrix->GetName(), matrix->GetTitle(), rows);
    info.g4OpticalProperties[matrix] = g4;
  }
  return g4;
}

static G4OpticalSurfaceFinish geant4_surface_finish(TGeoOpticalSurface::ESurfaceFinish f)   {
#define TO_G4_FINISH(x)  case TGeoOpticalSurface::kF##x : return x;
  switch(f)   {
    TO_G4_FINISH(polished);              // smooth perfectly polished surface
    TO_G4_FINISH(polishedfrontpainted);  // smooth top-layer (front) paint
    TO_G4_FINISH(polishedbackpainted);   // same is 'polished' but with a back-paint
 
    TO_G4_FINISH(ground);                // rough surface
    TO_G4_FINISH(groundfrontpainted);    // rough top-layer (front) paint
    TO_G4_FINISH(groundbackpainted);     // same as 'ground' but with a back-paint

    TO_G4_FINISH(polishedlumirrorair);   // mechanically polished surface, with lumirror
    TO_G4_FINISH(polishedlumirrorglue);  // mechanically polished surface, with lumirror & meltmount
    TO_G4_FINISH(polishedair);           // mechanically polished surface
    TO_G4_FINISH(polishedteflonair);     // mechanically polished surface, with teflon
    TO_G4_FINISH(polishedtioair);        // mechanically polished surface, with tio paint
    TO_G4_FINISH(polishedtyvekair);      // mechanically polished surface, with tyvek
    TO_G4_FINISH(polishedvm2000air);     // mechanically polished surface, with esr film
    TO_G4_FINISH(polishedvm2000glue);    // mechanically polished surface, with esr film & meltmount

    TO_G4_FINISH(etchedlumirrorair);     // chemically etched surface, with lumirror
    TO_G4_FINISH(etchedlumirrorglue);    // chemically etched surface, with lumirror & meltmount
    TO_G4_FINISH(etchedair);             // chemically etched surface
    TO_G4_FINISH(etchedteflonair);       // chemically etched surface, with teflon
    TO_G4_FINISH(etchedtioair);          // chemically etched surface, with tio paint
    TO_G4_FINISH(etchedtyvekair);        // chemically etched surface, with tyvek
    TO_G4_FINISH(etchedvm2000air);       // chemically etched surface, with esr film
    TO_G4_FINISH(etchedvm2000glue);      // chemically etched surface, with esr film & meltmount

    TO_G4_FINISH(groundlumirrorair);     // rough-cut surface, with lumirror
    TO_G4_FINISH(groundlumirrorglue);    // rough-cut surface, with lumirror & meltmount
    TO_G4_FINISH(groundair);             // rough-cut surface
    TO_G4_FINISH(groundteflonair);       // rough-cut surface, with teflon
    TO_G4_FINISH(groundtioair);          // rough-cut surface, with tio paint
    TO_G4_FINISH(groundtyvekair);        // rough-cut surface, with tyvek
    TO_G4_FINISH(groundvm2000air);       // rough-cut surface, with esr film
    TO_G4_FINISH(groundvm2000glue);      // rough-cut surface, with esr film & meltmount

    // for DAVIS model
    TO_G4_FINISH(Rough_LUT);             // rough surface
    TO_G4_FINISH(RoughTeflon_LUT);       // rough surface wrapped in Teflon tape
    TO_G4_FINISH(RoughESR_LUT);          // rough surface wrapped with ESR
    TO_G4_FINISH(RoughESRGrease_LUT);    // rough surface wrapped with ESR and coupled with opical grease
    TO_G4_FINISH(Polished_LUT);          // polished surface
    TO_G4_FINISH(PolishedTeflon_LUT);    // polished surface wrapped in Teflon tape
    TO_G4_FINISH(PolishedESR_LUT);       // polished surface wrapped with ESR
    TO_G4_FINISH(PolishedESRGrease_LUT); // polished surface wrapped with ESR and coupled with opical grease
    TO_G4_FINISH(Detector_LUT);          // polished surface with optical grease
  default:
    printout(ERROR,"Geant4Surfaces","++ Unknown finish style: %d [%s]. Assume polished!",
             int(f), TGeoOpticalSurface::FinishToString(f));
    return polished;
  }
#undef TO_G4_FINISH
}

static G4SurfaceType geant4_surface_type(TGeoOpticalSurface::ESurfaceType t)   {
#define TO_G4_TYPE(x)  case TGeoOpticalSurface::kT##x : return x;
  switch(t)   {
    TO_G4_TYPE(dielectric_metal);      // dielectric-metal interface
    TO_G4_TYPE(dielectric_dielectric); // dielectric-dielectric interface
    TO_G4_TYPE(dielectric_LUT);        // dielectric-Look-Up-Table interface
    TO_G4_TYPE(dielectric_LUTDAVIS);   // dielectric-Look-Up-Table DAVIS interface
    TO_G4_TYPE(dielectric_dichroic);   // dichroic filter interface
    TO_G4_TYPE(firsov);                // for Firsov Process
    TO_G4_TYPE(x_ray);                  // for x-ray mirror process
  default:
    printout(ERROR,"Geant4Surfaces","++ Unknown surface type: %d [%s]. Assume dielectric_metal!",
             int(t), TGeoOpticalSurface::TypeToString(t));
    return dielectric_metal;
  }
#undef TO_G4_TYPE
}

static G4OpticalSurfaceModel geant4_surface_model(TGeoOpticalSurface::ESurfaceModel surfMod)   {
#define TO_G4_MODEL(x)  case TGeoOpticalSurface::kM##x : return x;
  switch(surfMod)   {
    TO_G4_MODEL(glisur);   // original GEANT3 model
    TO_G4_MODEL(unified);  // UNIFIED model
    TO_G4_MODEL(LUT);      // Look-Up-Table model
    TO_G4_MODEL(DAVIS);    // DAVIS model
    TO_G4_MODEL(dichroic); // dichroic filter
  default:
    printout(ERROR,"Geant4Surfaces","++ Unknown surface model: %d [%s]. Assume glisur!",
             int(surfMod), TGeoOpticalSurface::ModelToString(surfMod));
    return glisur;
  }
#undef TO_G4_MODEL
}

/// Convert the optical surface to Geant4
void* Geant4Converter::handleOpticalSurface(TObject* surface) const    {
  TGeoOpticalSurface* optSurf    = (TGeoOpticalSurface*)surface;
  Geant4GeometryInfo& info = data();
  G4OpticalSurface*   g4   = info.g4OpticalSurfaces[optSurf];
  if ( !g4 ) {
    G4SurfaceType          type   = geant4_surface_type(optSurf->GetType());
    G4OpticalSurfaceModel  model  = geant4_surface_model(optSurf->GetModel());
    G4OpticalSurfaceFinish finish = geant4_surface_finish(optSurf->GetFinish());
    std::string name = make_NCName(optSurf->GetName());
    PrintLevel lvl = debugSurfaces ? ALWAYS : DEBUG;
    g4 = new G4OpticalSurface(name, model, finish, type, optSurf->GetValue());
    g4->SetSigmaAlpha(optSurf->GetSigmaAlpha());
    g4->SetPolish(optSurf->GetPolish());

    printout(lvl, "Geant4Converter",
             "++ Created OpticalSurface: %-18s type:%s model:%s finish:%s SigmaAlphs: %.3e Polish: %.3e",
             optSurf->GetName(),
             TGeoOpticalSurface::TypeToString(optSurf->GetType()),
             TGeoOpticalSurface::ModelToString(optSurf->GetModel()),
             TGeoOpticalSurface::FinishToString(optSurf->GetFinish()),
             optSurf->GetSigmaAlpha(), optSurf->GetPolish());
    ///
    /// Convert non-scalar properties from GDML tables
    G4MaterialPropertiesTable* tab = nullptr;
    TListIter itp(&optSurf->GetProperties());
    for(TObject* obj = itp.Next(); obj; obj = itp.Next())  {
      std::string exc_str;
      TNamed*      named  = (TNamed*)obj;
      TGDMLMatrix* matrix = info.manager->GetGDMLMatrix(named->GetTitle());
      const char*  cptr   = ::strstr(matrix->GetName(), GEANT4_TAG_IGNORE);
      if ( nullptr != cptr )  // Check if the property should not be passed to Geant4
        continue;

      if ( nullptr == tab )  {
        tab = new G4MaterialPropertiesTable();
        g4->SetMaterialPropertiesTable(tab);
      }

      Geant4GeometryInfo::PropertyVector* v =
        (Geant4GeometryInfo::PropertyVector*)handleMaterialProperties(matrix);
      if ( !v )  {  // Error!
        except("Geant4OpticalSurface","++ Failed to convert opt.surface %s. Property table %s is not defined!",
               optSurf->GetName(), named->GetTitle());
      }
      int idx = -1;
      try   {
        idx = tab->GetPropertyIndex(named->GetName());
      }
      catch(const std::exception& e)   {
        exc_str = e.what();
      }
      catch(...)   {
      }
      if ( idx < 0 )   {
        printout(ERROR, "Geant4Converter",
                 "++ UNKNOWN Geant4 Property: %-20s %s [IGNORED]",
                 exc_str.c_str(), named->GetName());
        continue;
      }
      // We need to convert the property from TGeo units to Geant4 units
      auto conv = g4PropertyConversion(idx);
      std::vector<double> bins(v->bins), vals(v->values);
      for(std::size_t i=0, count=v->bins.size(); i<count; ++i)
        bins[i] *= conv.first, vals[i] *= conv.second;
      G4MaterialPropertyVector* vec = new G4MaterialPropertyVector(&bins[0], &vals[0], bins.size());
      tab->AddProperty(named->GetName(), vec);
      
      printout(lvl, "Geant4Converter",
               "++       Property: %-20s [%ld x %ld] -->  %s",
               named->GetName(), matrix->GetRows(), matrix->GetCols(), named->GetTitle());
      for(std::size_t i=0, count=v->bins.size(); i<count; ++i)
        printout(lvl, named->GetName(),
                 "  Geant4: %8.3g [MeV]  TGeo: %8.3g [GeV] Conversion: %8.3g",
                 bins[i], v->bins[i], conv.first);
    }
    ///
    /// Convert scalar properties
#if ROOT_VERSION_CODE >= ROOT_VERSION(6,31,1)
    TListIter itc(&optSurf->GetConstProperties());
    for(TObject* obj = itc.Next(); obj; obj = itc.Next())  {
      std::string  exc_str;
      TNamed* named  = (TNamed*)obj;
      const char* cptr = ::strstr(named->GetName(), GEANT4_TAG_IGNORE);
      if ( nullptr != cptr )   {
        printout(INFO, name, "++ Ignore CONST property %s [%s].",
                 named->GetName(), named->GetTitle());
        continue;
      }
      cptr = ::strstr(named->GetTitle(), GEANT4_TAG_IGNORE);
      if ( nullptr != cptr )   {
        printout(INFO, name,"++ Ignore CONST property %s [%s].",
                 named->GetName(), named->GetTitle());
        continue;
      }
      Bool_t   err = kFALSE;
      Double_t value = info.manager->GetProperty(named->GetTitle(),&err);
      if ( err != kFALSE )   {
        except(name,
               "++ FAILED to create G4 material %s [Cannot convert const property: %s]",
               optSurf->GetName(), named->GetName());
      }
      if ( nullptr == tab )  {
        tab = new G4MaterialPropertiesTable();
        g4->SetMaterialPropertiesTable(tab);
      }
      int idx = -1;
      try   {
        idx = tab->GetConstPropertyIndex(named->GetName());
      }
      catch(const std::exception& e)   {
        exc_str = e.what();
      }
      catch(...)   {
      }
      if ( idx < 0 )   {
        printout(ERROR, name,
                 "++ UNKNOWN Geant4 CONST Property: %-20s %s [IGNORED]",
                 exc_str.c_str(), named->GetName());
        continue;
      }
      // We need to convert the property from TGeo units to Geant4 units
      double conv = g4ConstPropertyConversion(idx);
      printout(lvl, name, "++      CONST Property: %-20s %g * %g --> %g ",
               named->GetName(), value, conv, value * conv);
      tab->AddConstProperty(named->GetName(), value * conv);
    }
#endif  // ROOT_VERSION >= 6.31.1
    info.g4OpticalSurfaces[optSurf] = g4;
  }
  return g4;
}

/// Convert the skin surface to Geant4
void* Geant4Converter::handleSkinSurface(TObject* surface) const   {
  TGeoSkinSurface*    surf = (TGeoSkinSurface*)surface;
  Geant4GeometryInfo& info = data();
  G4LogicalSkinSurface* g4 = info.g4SkinSurfaces[surf];
  if ( !g4 ) {
    G4OpticalSurface* optSurf  = info.g4OpticalSurfaces[OpticalSurface(surf->GetSurface())];
    G4LogicalVolume*  v = info.g4Volumes[surf->GetVolume()];
    std::string name = make_NCName(surf->GetName());
    g4 = new G4LogicalSkinSurface(name, v, optSurf);
    printout(debugSurfaces ? ALWAYS : DEBUG, "Geant4Converter",
             "++ Created SkinSurface: %-18s  optical:%s",
             surf->GetName(), surf->GetSurface()->GetName());
    info.g4SkinSurfaces[surf] = g4;
  }
  return g4;
}

/// Convert the border surface to Geant4
void* Geant4Converter::handleBorderSurface(TObject* surface) const   {
  TGeoBorderSurface*    surf = (TGeoBorderSurface*)surface;
  Geant4GeometryInfo&   info = data();
  G4LogicalBorderSurface* g4 = info.g4BorderSurfaces[surf];
  if ( !g4 ) {
    G4OpticalSurface*  optSurf = info.g4OpticalSurfaces[OpticalSurface(surf->GetSurface())];
    G4VPhysicalVolume* n1 = info.g4Placements[surf->GetNode1()];
    G4VPhysicalVolume* n2 = info.g4Placements[surf->GetNode2()];
    std::string name = make_NCName(surf->GetName());
    g4 = new G4LogicalBorderSurface(name, n1, n2, optSurf);
    printout(debugSurfaces ? ALWAYS : DEBUG, "Geant4Converter",
             "++ Created BorderSurface: %-18s  optical:%s",
             surf->GetName(), surf->GetSurface()->GetName());
    info.g4BorderSurfaces[surf] = g4;
  }
  return g4;
}

/// Convert the geometry type SensitiveDetector into the corresponding Geant4 object(s).
void Geant4Converter::printSensitive(SensitiveDetector sens_det, const std::set<const TGeoVolume*>& /* volumes */) const {
  Geant4GeometryInfo&          info = data();
  std::set<const TGeoVolume*>& volset = info.sensitives[sens_det];
  SensitiveDetector            sd = sens_det;
  std::stringstream str;

  printout(INFO, "Geant4Converter", "++ SensitiveDetector: %-18s %-20s Hits:%-16s", sd.name(), ("[" + sd.type() + "]").c_str(),
           sd.hitsCollection().c_str());
  str << "                    | " << "Cutoff:" << std::setw(6) << std::left
      << sd.energyCutoff() << std::setw(5) << std::right << volset.size()
      << " volumes ";
  if (sd.region().isValid())
    str << " Region:" << std::setw(12) << std::left << sd.region().name();
  if (sd.limits().isValid())
    str << " Limits:" << std::setw(12) << std::left << sd.limits().name();
  str << ".";
  printout(INFO, "Geant4Converter", str.str().c_str());

  for (const auto i : volset )  {
    std::map<Volume, G4LogicalVolume*>::iterator v = info.g4Volumes.find(i);
    if ( v != info.g4Volumes.end() )   {
      G4LogicalVolume* vol = (*v).second;
      str.str("");
      str << "                                   | " << "Volume:" << std::setw(24) << std::left << vol->GetName() << " "
          << vol->GetNoDaughters() << " daughters.";
      printout(INFO, "Geant4Converter", str.str().c_str());
    }
  }
}

std::string printSolid(G4VSolid* sol) {
  std::stringstream str;
  if (typeid(*sol) == typeid(G4Box)) {
    const G4Box* b = (G4Box*) sol;
    str << "++ Box: x=" << b->GetXHalfLength() << " y=" << b->GetYHalfLength() << " z=" << b->GetZHalfLength();
  }
  else if (typeid(*sol) == typeid(G4Tubs)) {
    const G4Tubs* t = (const G4Tubs*) sol;
    str << " Tubs: Ri=" << t->GetInnerRadius() << " Ra=" << t->GetOuterRadius() << " z/2=" << t->GetZHalfLength() << " Phi="
        << t->GetStartPhiAngle() << "..." << t->GetDeltaPhiAngle();
  }
  return str.str();
}

/// Print G4 placement
void* Geant4Converter::printPlacement(const std::string& name, const TGeoNode* node) const {
  Geant4GeometryInfo& info = data();
  G4VPhysicalVolume*  g4   = info.g4Placements[node];
  G4LogicalVolume*    vol  = info.g4Volumes[node->GetVolume()];
  G4LogicalVolume*    mot  = info.g4Volumes[node->GetMotherVolume()];
  G4VSolid*           sol  = vol->GetSolid();
  G4ThreeVector       tr   = g4->GetObjectTranslation();
  G4VSensitiveDetector* sd = vol->GetSensitiveDetector();
  if ( !sd )  {
    return g4;
  }
  std::stringstream str;
  str << "G4Cnv::placement: + " << name << " No:" << node->GetNumber() << " Vol:" << vol->GetName() << " Solid:"
      << sol->GetName();
  printout(outputLevel, "G4Placement", str.str().c_str());
  str.str("");
  str << "                  |" << " Loc: x=" << tr.x() << " y=" << tr.y() << " z=" << tr.z();
  printout(outputLevel, "G4Placement", str.str().c_str());
  printout(outputLevel, "G4Placement", printSolid(sol).c_str());
  str.str("");
  str << "                  |" << " Ndau:" << vol->GetNoDaughters()
      << " physvols." << " Mat:" << vol->GetMaterial()->GetName()
      << " Mother:" << (char*) (mot ? mot->GetName().c_str() : "---");
  printout(outputLevel, "G4Placement", str.str().c_str());
  str.str("");
  str << "                  |" << " SD:" << sd->GetName();
  printout(outputLevel, "G4Placement", str.str().c_str());
  return g4;
}

/// Create geometry conversion
Geant4Converter& Geant4Converter::create(DetElement top) {
  typedef std::map<const TGeoNode*, std::vector<TGeoNode*> > _DAU;
  TTimeStamp start;
  _DAU daughters;
  Geant4GeometryInfo& geo = this->init();
  World wrld = top.world();

  m_data->clear();
  m_set_data->clear();
  m_daughters = &daughters;
  geo.manager = &wrld.detectorDescription().manager();
  this->collect(top, geo);
  this->checkOverlaps = false;
  // We do not have to handle defines etc.
  // All positions and the like are not really named.
  // Hence, start creating the G4 objects for materials, solids and log volumes.
  handleArray(this, geo.manager->GetListOfGDMLMatrices(), &Geant4Converter::handleMaterialProperties);
  handleArray(this, geo.manager->GetListOfOpticalSurfaces(), &Geant4Converter::handleOpticalSurface);
  
  handle(this,     geo.volumes, &Geant4Converter::collectVolume);
  handle(this,     geo.solids,  &Geant4Converter::handleSolid);
  printout(outputLevel, "Geant4Converter", "++ Handled %ld solids.", geo.solids.size());
  handleRefs(this, geo.vis,     &Geant4Converter::handleVis);
  printout(outputLevel, "Geant4Converter", "++ Handled %ld visualization attributes.", geo.vis.size());
  handleMap(this,  geo.limits,  &Geant4Converter::handleLimitSet);
  printout(outputLevel, "Geant4Converter", "++ Handled %ld limit sets.", geo.limits.size());
  handleMap(this,  geo.regions, &Geant4Converter::handleRegion);
  printout(outputLevel, "Geant4Converter", "++ Handled %ld regions.", geo.regions.size());
  handle(this,     geo.volumes, &Geant4Converter::handleVolume);
  printout(outputLevel, "Geant4Converter", "++ Handled %ld volumes.", geo.volumes.size());
  handleRMap(this, *m_data,     &Geant4Converter::handleAssembly);
  // Now place all this stuff appropriately
  //handleRMap(this, *m_data,     &Geant4Converter::handlePlacement);
  std::map<int, std::vector<const TGeoNode*> >::const_reverse_iterator i = m_data->rbegin();
  for ( ; i != m_data->rend(); ++i )  {
    for ( const TGeoNode* node : i->second )  {
      this->handlePlacement(node->GetName(), node);
    }
  }
  /// Handle concrete surfaces
  handleArray(this, geo.manager->GetListOfSkinSurfaces(),   &Geant4Converter::handleSkinSurface);
  handleArray(this, geo.manager->GetListOfBorderSurfaces(), &Geant4Converter::handleBorderSurface);
  //==================== Fields
  handleProperties(m_detDesc.properties());
  if ( printSensitives )  {
    handleMap(this, geo.sensitives, &Geant4Converter::printSensitive);
  }
  if ( printPlacements )  {
    handleRMap(this, *m_data, &Geant4Converter::printPlacement);
  }

  m_daughters = nullptr;
  geo.setWorld(top.placement().ptr());
  geo.valid = true;
  TTimeStamp stop;
  printout(INFO, "Geant4Converter",
           "+++  Successfully converted geometry to Geant4. [%7.3f seconds]",
           stop.AsDouble()-start.AsDouble() );
  return *this;
}