CADGeometryUtils.cxx

Go to the documentation of this file.

// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
// All rights not expressly granted are reserved.
//
// This software is distributed under the terms of the GNU General Public
// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
//
// In applying this license CERN does not waive the privileges and immunities
// granted to it by virtue of its status as an Intergovernmental Organization
// or submit itself to any jurisdiction.
 
#include "CADSupport/CADGeometryUtils.h"
#include "DetectorsBase/MaterialManager.h"
#include <CommonUtils/FileSystemUtils.h>
#include <TGeoVolume.h>
#include <TGeoNode.h>
#include <TGeoMaterial.h>
#include <TGeoMedium.h>
#include <TInterpreter.h>
#include <TROOT.h>
#include <TString.h>
#include <TGlobal.h>
#include <fairlogger/Logger.h>
#include <atomic>
#include <cctype>
#include <filesystem>
#include <fstream>
#include <functional>
#include <string>
#include <unordered_map>
#include <unordered_set>
#include <vector>
 
namespace o2::cad
{
 
TGeoVolume* buildCADVolumeFromMacro(const std::string& macroFile, const std::string& instanceTag)
{
  if (macroFile.empty()) {
    return nullptr;
  }
  auto expandedHookFileName = o2::utils::expandShellVarsInFileName(macroFile);
  if (!std::filesystem::exists(expandedHookFileName)) {
    LOG(error) << "External geometry macro " << expandedHookFileName << " does not exist";
    return nullptr;
  }
 
  // JIT each macro into its own namespace, since every converter macro defines the same symbols; includes stay global.
  std::ifstream macroStream(expandedHookFileName, std::ios::in);
  if (!macroStream.is_open()) {
    LOG(error) << "Cannot open external geometry macro " << expandedHookFileName;
    return nullptr;
  }
  std::string preamble; // #include (and other top-level preprocessor) lines -> global scope
  std::string body;     // everything else -> wrapped into a unique namespace
  std::string line;
  while (std::getline(macroStream, line)) {
    auto firstNonSpace = line.find_first_not_of(" \t");
    if (firstNonSpace != std::string::npos && line[firstNonSpace] == '#') {
      preamble += line + "\n";
    } else {
      body += line + "\n";
    }
  }
 
  // build a unique, valid C++ identifier for the namespace
  static std::atomic<int> instanceCounter{0};
  std::string ns = std::string("o2_cadgeom_") + instanceTag + "_" + std::to_string(instanceCounter++);
  for (auto& c : ns) {
    if (!std::isalnum(static_cast<unsigned char>(c)) && c != '_') {
      c = '_';
    }
  }
 
  const std::string wrapped = preamble + "\nnamespace " + ns + " {\n" + body + "\n}\n";
  if (!gInterpreter->Declare(wrapped.c_str())) {
    LOG(error) << "Failed to JIT external geometry macro " << expandedHookFileName;
    return nullptr;
  }
 
  // retrieve the builder hook from the unique namespace
  const std::string globalName = "__" + ns + "_hook__";
  gROOT->ProcessLine(Form("std::function<TGeoVolume*()> %s = %s::get_builder_hook_unchecked();",
                          globalName.c_str(), ns.c_str()));
  auto global = gROOT->GetGlobal(globalName.c_str());
  if (!global) {
    LOG(error) << "Could not retrieve geometry builder hook from macro " << expandedHookFileName;
    return nullptr;
  }
  auto hook = *reinterpret_cast<std::function<TGeoVolume*()>*>(global->GetAddress());
  LOG(info) << "CAD geometry hook initialized from file " << expandedHookFileName << " (namespace " << ns << ")";
 
  auto top = hook();
  if (!top) {
    LOG(error) << "CAD geometry macro " << expandedHookFileName << " did not return a top volume";
  }
  return top;
}
 
void remapCADMedia(TGeoVolume* top, const char* modulename)
{
  std::unordered_map<TGeoMedium*, TGeoMedium*> medium_ptr_mapping;
  std::unordered_set<TGeoVolume*> volumes_already_treated;
  // a material may back several media (the `_NF` twins), so materials are deduplicated apart from media
  std::unordered_map<std::string, int> material_index;
  int counter = 1;
  int matcounter = 1;
 
  // The transformer function
  auto transform_media = [&](TGeoVolume* vol_) {
    if (volumes_already_treated.find(vol_) != volumes_already_treated.end()) {
      // this volume was already transformed
      return;
    }
    volumes_already_treated.insert(vol_);
 
    if (dynamic_cast<TGeoVolumeAssembly*>(vol_)) {
      // do nothing for assemblies (they don't have a medium)
      return;
    }
 
    auto medium = vol_->GetMedium();
    if (!medium) {
      return;
    }
 
    auto iter = medium_ptr_mapping.find(medium);
    if (iter != medium_ptr_mapping.end()) {
      // This medium has already been transformed, so
      // we just update the volume
      vol_->SetMedium(iter->second);
      return;
    } else {
      LOG(info) << "Transforming media with name " << medium->GetName() << " for volume " << vol_->GetName();
 
      // we found a medium, not yet treated
      auto curr_mat = medium->GetMaterial();
      auto& matmgr = o2::base::MaterialManager::Instance();
 
      // Register the material once, however many media wear it.
      const std::string matname(curr_mat->GetName());
      auto itmat = material_index.find(matname);
      int imat;
      if (itmat != material_index.end()) {
        imat = itmat->second;
      } else {
        imat = matcounter++;
        // A TGeoMixture goes through Mixture() so Geant keeps its element composition.
        if (auto* mix = dynamic_cast<TGeoMixture*>(curr_mat)) {
          const Int_t nel = mix->GetNelements();
          std::vector<Float_t> a(nel), z(nel), w(nel);
          for (Int_t i = 0; i < nel; ++i) {
            a[i] = mix->GetAmixt()[i];
            z[i] = mix->GetZmixt()[i];
            w[i] = mix->GetWmixt()[i];
          }
          matmgr.Mixture(modulename, imat, curr_mat->GetName(), a.data(), z.data(),
                         curr_mat->GetDensity(), nel, w.data());
        } else {
          matmgr.Material(modulename, imat, curr_mat->GetName(), curr_mat->GetA(), curr_mat->GetZ(), curr_mat->GetDensity(), curr_mat->GetRadLen(), curr_mat->GetIntLen());
        }
        material_index[matname] = imat;
      }
      // TGeo medium params are stored in a flat array with the following convention
      // fParams[0] = isvol;
      // fParams[1] = ifield;
      // fParams[2] = fieldm;
      // fParams[3] = tmaxfd;
      // fParams[4] = stemax;
      // fParams[5] = deemax;
      // fParams[6] = epsil;
      // fParams[7] = stmin;
      const auto isvol = medium->GetParam(0);
      const auto isxfld = medium->GetParam(1);
      const auto sxmgmx = medium->GetParam(2);
      const auto tmaxfd = medium->GetParam(3);
      const auto stemax = medium->GetParam(4);
      const auto deemax = medium->GetParam(5);
      const auto epsil = medium->GetParam(6);
      const auto stmin = medium->GetParam(7);
 
      matmgr.Medium(modulename, counter, medium->GetName(), imat, isvol, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
 
      // there will be new Material and Medium objects; fetch them
      auto new_med = matmgr.getTGeoMedium(modulename, counter);
 
      // insert into cache
      medium_ptr_mapping[medium] = new_med;
      vol_->SetMedium(new_med);
      counter++;
    }
  }; // end transformer lambda
 
  // a generic volume walker
  std::function<void(TGeoVolume*)> visit_volume;
  visit_volume = [&](TGeoVolume* vol) -> void {
    if (!vol) {
      return;
    }
 
    // call the transformer
    transform_media(vol);
 
    // Recurse into daughters
    const int nd = vol->GetNdaughters();
    for (int i = 0; i < nd; ++i) {
      TGeoNode* node = vol->GetNode(i);
      if (!node) {
        continue;
      }
      TGeoVolume* child = node->GetVolume();
      if (!child) {
        continue;
      }
 
      visit_volume(child);
    }
  };
 
  visit_volume(top);
}
 
} // namespace o2::cad