KfSetupBuilder.cxx

Go to the documentation of this file.

/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
   SPDX-License-Identifier: GPL-3.0-only
   Authors: Sergei Zharko [committer] */

 
#include "KfSetupBuilder.h"
 
#include <boost/filesystem.hpp>
 
#include <numeric>
 
using cbm::algo::kf::SetupBuilder;
 
namespace fs = boost::filesystem;
 
// ---------------------------------------------------------------------------------------------------------------------
//
std::string SetupBuilder::LayerProperty::ToString() const
{
  return fmt::format(
    "GeoLayer: detId={}, locId={}, zRef={:9.4}cm, z=[{:9.4}, {:9.4}]cm, xFull=[{:9.4}, {:9.4}]cm, "
    "yFull=[{:9.4, {:9.4}}]cm, xAct=[{:9.4}, {:9.4}]cm, yAct=[{:9.4}, {:9.4}]cm, fieldType={}, timeInfo={}",
    detId, locId, zRef, zMin, zMax, xMinFull, xMaxFull, yMinFull, yMaxFull, xMinAct, xMaxAct, yMinAct, yMaxAct,
    static_cast<int>(fieldType), timeInfo);
}
 
// ---------------------------------------------------------------------------------------------------------------------
//
void SetupBuilder::AddLayer(SetupBuilder::LayerProperty geoLayer)
{
  fbReady                     = false;
  auto [itAdded, wasInserted] = fGeoLayers.insert(geoLayer);
  if (!wasInserted) {
    throw std::logic_error(fmt::format(
      "SetupBuilder::AddLayer: attempt of adding a duplicating geometry layer: {}\n"
      "The next layers were already added: {}",
      itAdded->ToString(),
      std::accumulate(fGeoLayers.begin(), fGeoLayers.end(), std::string{}, [](std::string s, const auto& l) {
        return std::move(s) + "\n\t - " + std::move(l.ToString());
      })));
  }
}
 
 
// ---------------------------------------------------------------------------------------------------------------------
//
void SetupBuilder::Init()
{
  LOG(info) << "kf::SetupBuilder initialization: ...";
 
  if (!fbIfFieldFunctionSet) {
    throw std::runtime_error("kf::SetupBuilder: no field is provided");
  }
  if (fbIfSetFromSetup) {
    LOG(info) << "kf::SetupBuilder initialization: done";
    fbReady = true;
    return;
  }
  else {
    if (fGeoLayers.size() == 0) {
      throw std::runtime_error("kf::SetupBuilder: no geometry layers initialized");
    }
    if (!fpMaterialMapFactory.get()) {
      throw std::runtime_error("kf::SetupBuilder: no material map factory provided");
    }
    if (!fbIfTargetSet) {
      throw std::runtime_error("kf::SetupBuilder: target properties were not set");
    }
  }
 
  bool bMaterialLoaded{false};
  if (!fsMaterialCacheFile.empty()) {
    bMaterialLoaded = this->LoadMaterial();
  }
 
  // Add material to the target:
  if (!bMaterialLoaded) {
    // Material budget
    double zMin{fTarget.GetZ() + kTargetCenterOffset};
    double zMax{fTarget.GetZ() + fTarget.GetDz() * kTargetMaterialOffset};
    double rMax{fTarget.GetR() * kTargetMaterialTransverseSizeMargin};
    fTarget.SetMaterial(
      fpMaterialMapFactory->GenerateMaterialMap(fTarget.GetZ(), zMin, zMax, rMax, kTargetMaterialMapNofBins));
  }
 
  // Init geometry layers
  {
    // Estimate acceptance slope
    double acceptanceSlope = -1.;  // arbitrary negative value
    for (const auto& layer : fGeoLayers) {
      double xyMax = std::max(std::max(std::fabs(layer.xMaxFull), std::fabs(layer.xMinFull)),
                              std::max(std::fabs(layer.yMaxFull), std::fabs(layer.yMinFull)));
      double tCurr = xyMax / (layer.zRef - fTarget.GetZ());
      assert(tCurr > 0);
      acceptanceSlope = std::max(acceptanceSlope, tCurr);
    }
 
    if (!bMaterialLoaded) {
      fvMaterial.clear();
      fvMaterial.reserve(fGeoLayers.size());
    }
    double zLast{fTarget.GetZ() + fTarget.GetDz() * kTargetMaterialOffset};
    for (auto layerIt = fGeoLayers.cbegin(); layerIt != fGeoLayers.cend(); ++layerIt) {
      double z1 = layerIt->zMax;
      double z2 = z1;
      if (std::next(layerIt) != fGeoLayers.cend()) {
        z2 = std::next(layerIt)->zMin;
      }
      double zNew{0.5 * (z1 + z2)};
      double xyMax{acceptanceSlope * (layerIt->zRef - fTarget.GetZ())};
      if (!bMaterialLoaded) {
        auto material{fpMaterialMapFactory->GenerateMaterialMap(layerIt->zRef, zLast, zNew, xyMax, fMatMapNofBins)};
        fvMaterial.push_back(std::move(material));
      }
 
      // Note: square material maps to follow the material budget map regions
      // TODO: test field xMax, xMin
      double fieldXmax = std::max(std::fabs(layerIt->xMaxFull), std::fabs(layerIt->xMinFull));
      double fieldYmax = std::max(std::fabs(layerIt->yMaxFull), std::fabs(layerIt->yMinFull));
      fFieldBuilder.AddSliceReference(fieldXmax, fieldYmax, layerIt->zRef);
      fModuleIndexFactory.AddComponent<int>(layerIt->detId, layerIt->locId, layerIt->zRef);
      zLast = zNew;
    }
  }
 
  if (!bMaterialLoaded && !fsMaterialCacheFile.empty()) {
    this->StoreMaterial();
  }
 
  fbIfGeoLayersInit = true;
  fbReady           = true;
  LOG(info) << "kf::SetupBuilder initialization: done";
}
 
// ---------------------------------------------------------------------------------------------------------------------
//
std::string SetupBuilder::InitStatusMsg() const
{
  return fmt::format(                           //
    "kf::SetupBuilder initialization: failed:"  //
    "\n - target property set:      {}"         //
    "\n - geo layers added:         {}"         //
    "\n - field function set:       {}"         //
    "\n - set from setup:           {}"         //
    "\n - material map creator set: {}",        //
    fbIfTargetSet, fbIfGeoLayersInit, fbIfFieldFunctionSet, fbIfSetFromSetup, nullptr != fpMaterialMapFactory);
}
 
// ---------------------------------------------------------------------------------------------------------------------
//
bool SetupBuilder::LoadMaterial()
{
  auto path = fs::absolute(fs::weakly_canonical(fsMaterialCacheFile));
 
  std::ifstream ifs(path.string(), std::ios::binary);
  if (!ifs) {  // File does not exist yet
    return false;
  }
  try {
    boost::archive::binary_iarchive ia(ifs);
    MaterialMap targetMat;
    size_t refHash;
    ia >> refHash;
    if (refHash != fGeoHash) {
      LOG(warn) << "kf::SetupBuilder::LoadMaterial: reference hash from input file \"" << path.string()
                << "\" "
                   "diverges from one, obtained from the actual detector setup geometry. Material budget will be "
                   "re-generated, and a new file will be created";
      return false;
    }
    ia >> targetMat;
    ia >> fvMaterial;
    fTarget.SetMaterial(targetMat);
  }
  catch (const std::exception& err) {
    LOG(warn) << "kf::SetupBuilder::LoadMaterial: input file \"" << path.string()
              << "\" has inconsistent format "
                 "or was corrupted. The material maps will be re-generated";
    return false;
  }
  LOG(info) << "kf::SetupBuilder::LoadMaterial: the material maps were loaded from cache file \"" << path.string()
            << "\"";
  return true;
}
 
// ---------------------------------------------------------------------------------------------------------------------
//
void SetupBuilder::Reset()
{
  fbReady              = false;
  fbIfTargetSet        = false;
  fbIfGeoLayersInit    = false;
  fbIfFieldFunctionSet = false;
  fbIfSetFromSetup     = false;
  fFieldBuilder.Reset();
  fModuleIndexFactory.Reset();
  fvMaterial.clear();
  fGeoLayers.clear();
  fpMaterialMapFactory = nullptr;
}
 
// ---------------------------------------------------------------------------------------------------------------------
//
void SetupBuilder::SetTargetProperty(double x, double y, double z, double dz, double r)
{
  fbReady = false;
  fTarget.SetX(x);
  fTarget.SetY(y);
  fTarget.SetZ(z);
  fTarget.SetDz(dz);
  fTarget.SetR(r);
  fFieldBuilder.SetTarget(x, y, z);
  fbIfTargetSet = true;
}
 
// ---------------------------------------------------------------------------------------------------------------------
//
void SetupBuilder::Store(const Setup<double>& setup, const std::string& fileName)
{
  auto path = fs::absolute(fs::weakly_canonical(fileName));
  fs::create_directories(path.parent_path());
 
  std::ofstream ofs(path.string(), std::ios::binary);
  if (!ofs) {
    throw std::runtime_error(
      fmt::format("kf::SetupBuilder::Store: failed opening file \"{}\" to store the setup", path.string()));
  }
  LOG(info) << "kf::SetupBuilder: writing setup into file \"" << path.string() << '\"';
  boost::archive::binary_oarchive oa(ofs);
  oa << setup;
}
 
// ---------------------------------------------------------------------------------------------------------------------
//
void SetupBuilder::StoreMaterial() const
{
  auto path = fs::absolute(fs::weakly_canonical(fsMaterialCacheFile));
  fs::create_directories(path.parent_path());
 
  std::ofstream ofs(path.string(), std::ios::binary);
  if (!ofs) {
    throw std::runtime_error(
      fmt::format("kf::SetupBuilder::Store: failed opening file \"{}\" to store the material cache", path.string()));
  }
  boost::archive::binary_oarchive oa(ofs);
  oa << fGeoHash;
  oa << fTarget.GetMaterial();
  oa << fvMaterial;
}