Geant4Converter.cpp

    }
    data().g4Limits[limitset] = g4;
  }
  return g4;
}

/// Convert the geometry visualisation attributes to the corresponding Geant4 object(s).
void* Geant4Converter::handleVis(const 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 {
  map < string, string > processors;
  static int s_idd = 9999999;
  string id;
  for (Detector::Properties::const_iterator i = prp.begin(); i != prp.end(); ++i) {
    const string& nam = (*i).first;
    const Detector::PropertyValues& vals = (*i).second;
    if (nam.substr(0, 6) == "geant4") {
      Detector::PropertyValues::const_iterator id_it = vals.find("id");
      if (id_it != vals.end()) {
        id = (*id_it).second;
      }
      else {
        char txt[32];
        ::snprintf(txt, sizeof(txt), "%d", ++s_idd);
        id = txt;
      }
      processors.emplace(id, nam);
    }
  }
  for (map<string, string>::const_iterator i = processors.begin(); i != processors.end(); ++i) {
    const GeoHandler* hdlr = this;
    string nam = (*i).second;
    const Detector::PropertyValues& vals = prp[nam];
    string type = vals.find("type")->second;
    string tag = type + "_Geant4_action";
    Detector* detPtr = const_cast<Detector*>(&m_detDesc);
    long result = PluginService::Create<long>(tag, detPtr, hdlr, &vals);
    if (0 == result) {
      throw runtime_error("Failed to locate plugin to interprete files of type"
                          " \"" + tag + "\" - no factory:" + type);
    }
    result = *(long*) result;
    if (result != 1) {
      throw runtime_error("Failed to invoke the plugin " + tag + " of type " + type);
    }
    printout(outputLevel, "Geant4Converter", "+++++ Executed Successfully Geant4 setup module *%s*.", type.c_str());
  }
}

#if ROOT_VERSION_CODE >= ROOT_VERSION(6,17,0)
/// Convert the geometry type material into the corresponding Geant4 object(s).
void* Geant4Converter::handleMaterialProperties(TObject* matrix) const    {
  TGDMLMatrix* gdmlMat = (TGDMLMatrix*)matrix;
  Geant4GeometryInfo& info = data();
  Geant4GeometryInfo::PropertyVector* g4 = info.g4OpticalProperties[gdmlMat];
  if (!g4) {
    PrintLevel lvl = debugMaterials ? ALWAYS : outputLevel;
    g4 = new Geant4GeometryInfo::PropertyVector();
    size_t rows = gdmlMat->GetRows();
    g4->name    = gdmlMat->GetName();
    g4->title   = gdmlMat->GetTitle();
    g4->bins.reserve(rows);
    g4->values.reserve(rows);
    for(size_t i=0; i<rows; ++i)  {
      g4->bins.emplace_back(gdmlMat->Get(i,0)  /*   *CLHEP::eV/units::eV   */);
      g4->values.emplace_back(gdmlMat->Get(i,1));
    }
    printout(lvl, "Geant4Converter", "++ Successfully converted material property:%s : %s [%ld rows]",
             gdmlMat->GetName(), gdmlMat->GetTitle(), rows);
    info.g4OpticalProperties[gdmlMat] = 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());
    g4 = new G4OpticalSurface(optSurf->GetName(), model, finish, type, optSurf->GetValue());
    g4->SetSigmaAlpha(optSurf->GetSigmaAlpha());
    // not implemented: g4->SetPolish(s->GetPolish());
    printout(debugSurfaces ? ALWAYS : DEBUG, "Geant4Converter",
             "++ Created OpticalSurface: %-18s type:%s model:%s finish:%s",
             optSurf->GetName(),
             TGeoOpticalSurface::TypeToString(optSurf->GetType()),
             TGeoOpticalSurface::ModelToString(optSurf->GetModel()),
             TGeoOpticalSurface::FinishToString(optSurf->GetFinish()));
    G4MaterialPropertiesTable* tab = 0;
    TListIter it(&optSurf->GetProperties());
    for(TObject* obj = it.Next(); obj; obj = it.Next())  {
      TNamed*      named = (TNamed*)obj;
      TGDMLMatrix* matrix = info.manager->GetGDMLMatrix(named->GetTitle());
      if ( 0 == 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 = tab->GetPropertyIndex(named->GetName(), false);
      if ( idx < 0 )   {
        printout(ERROR, "Geant4Converter", "++ UNKNOWN Geant4 Property: %-20s [IGNORED]", named->GetName());
        continue;
      }
      // We need to convert the property from TGeo units to Geant4 units
      auto conv = g4PropertyConversion(idx);
      vector<double> bins(v->bins), vals(v->values);
      for(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(debugSurfaces ? ALWAYS : DEBUG, "Geant4Converter",
               "++       Property: %-20s [%ld x %ld] -->  %s",
               named->GetName(), matrix->GetRows(), matrix->GetCols(), named->GetTitle());
      for(size_t i=0, count=v->bins.size(); i<count; ++i)
        printout(debugSurfaces ? ALWAYS : DEBUG, named->GetName(),
                 "  Geant4: %8.3g [MeV]  TGeo: %8.3g [GeV] Conversion: %8.3g",
                 bins[i], v->bins[i], conv.first);
    }
    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()];
    g4 = new G4LogicalSkinSurface(surf->GetName(), 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()];
    g4 = new G4LogicalBorderSurface(surf->GetName(), n1, n2, optSurf);
    printout(debugSurfaces ? ALWAYS : DEBUG, "Geant4Converter",
             "++ Created BorderSurface: %-18s  optical:%s",
             surf->GetName(), surf->GetSurface()->GetName());
    info.g4BorderSurfaces[surf] = g4;
  }
  return g4;
}
#endif

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

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

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

string printSolid(G4VSolid* sol) {
  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 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;
  }
  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;
}

template <typename O, typename C, typename F> void handleRefs(const O* o, const C& c, F pmf) {
  for (typename C::const_iterator i = c.begin(); i != c.end(); ++i) {
    //(o->*pmf)((*i)->GetName(), *i);
    (o->*pmf)("", *i);
  }
}

template <typename O, typename C, typename F> void handle(const O* o, const C& c, F pmf) {
  for (typename C::const_iterator i = c.begin(); i != c.end(); ++i) {
    (o->*pmf)((*i)->GetName(), *i);
  }
}

template <typename O, typename F> void handleArray(const O* o, const TObjArray* c, F pmf) {
  TObjArrayIter arr(c);
  for(TObject* i = arr.Next(); i; i=arr.Next())
    (o->*pmf)(i);
}

template <typename O, typename C, typename F> void handleMap(const O* o, const C& c, F pmf) {
  for (typename C::const_iterator i = c.begin(); i != c.end(); ++i)
    (o->*pmf)((*i).first, (*i).second);
}

template <typename O, typename C, typename F> void handleRMap(const O* o, const C& c, F pmf) {
  for (typename C::const_reverse_iterator i = c.rbegin(); i != c.rend(); ++i)  {
    //cout << "Handle RMAP [ " << (*i).first << " ]" << endl;
    handle(o, (*i).second, pmf);
  }
}

/// Create geometry conversion
Geant4Converter& Geant4Converter::create(DetElement top) {
  Geant4GeometryInfo& geo = this->init();
  World wrld = top.world();
  m_data->clear();
  geo.manager = &wrld.detectorDescription().manager();
  collect(top, geo);
  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.

  //outputLevel = WARNING;
  //setPrintLevel(VERBOSE);
  debugMaterials  = true;
  debugElements   = true;
  
#if ROOT_VERSION_CODE >= ROOT_VERSION(6,17,0)
  handleArray(this, geo.manager->GetListOfGDMLMatrices(), &Geant4Converter::handleMaterialProperties);
  handleArray(this, geo.manager->GetListOfOpticalSurfaces(), &Geant4Converter::handleOpticalSurface);
#endif
  
  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);
#if ROOT_VERSION_CODE >= ROOT_VERSION(6,17,0)
  /// Handle concrete surfaces
  handleArray(this, geo.manager->GetListOfSkinSurfaces(),   &Geant4Converter::handleSkinSurface);
  handleArray(this, geo.manager->GetListOfBorderSurfaces(), &Geant4Converter::handleBorderSurface);
#endif
  //==================== Fields
  handleProperties(m_detDesc.properties());
  if ( printSensitives )  {
    handleMap(this, geo.sensitives, &Geant4Converter::printSensitive);
  }
  if ( printPlacements )  {
    handleRMap(this, *m_data, &Geant4Converter::printPlacement);
  }

  geo.setWorld(top.placement().ptr());
  geo.valid = true;
  printout(INFO, "Geant4Converter", "+++  Successfully converted geometry to Geant4.");
  return *this;
}