CbmKfTrackingGeoSetupFactory.cxx

Go to the documentation of this file.

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

 
#include "CbmKfTrackingGeoSetupFactory.h"
 
#include "CbmKfTarget.h"
#include "CbmRecoSetupManager.h"
#include "CbmSetup.h"
#include "CbmStsFindTracks.h"
#include "KfDefs.h"
#include "KfMaterialMapFactory.h"
#include "KfToolsField.h"
 
#include <FairField.h>
#include <FairRootManager.h>
#include <FairRunAna.h>
#include <Logger.h>
 
#include <boost/filesystem.hpp>
 
#include <functional>
 
using cbm::algo::ca::EDetectorID;
using cbm::kf::TrackingGeoSetupFactory;
using kf::tools::MaterialMapFactory;
 
// FIXME: provide a validation function to ensure, that all necessary input is properly initialized:
// 1) Det. interfaces
// 2) Magnetic field
// 3) Target
// 4) GeoManager
 
// ---------------------------------------------------------------------------------------------------------------------
//
TrackingGeoSetupFactory::TrackingGeoSetupFactory(bool requireHitPresence)
try {
  using cbm::algo::kf::EFieldType;
  using cbm::algo::kf::GlobalField;
  using cbm::algo::kf::defs::MinField;
  LOG(info) << "cbm::kf::TrackingGeoSetupFactory: initializing ...";
 
  CheckDetectorPresence(requireHitPresence);
 
  // Reco setup instance
  if (!RecoSetupManager::Instance()->IsInitialized()) {
    throw std::runtime_error("cbm::RecoSetupManager not initialized");
  }
 
  // Magnetic field initialization
  cbm::kf::tools::SetOriginalCbmField();
 
  fBuilder.SetFieldFunction(GlobalField::fgOriginalField, GlobalField::fgOriginalFieldType);
 
  // NOTE: assigning a field type for a station: a field on the station is Normal, the field value in the middle
  //       of the station is above the threshold.
  auto GuessFieldType = [&](double z) -> EFieldType {
    auto [Bx, By, Bz] = GlobalField::fgOriginalField(0., 0., z);
    double absBsq     = Bx * Bx + By * By + Bz * Bz;
    return absBsq < MinField<double> * MinField<double> ? EFieldType::Null : EFieldType::Normal;
  };
 
  // Set tracking station properties
  auto CollectStations = [&](const auto* pIfs, EDetectorID detID) -> void {
    if (!fvbDetInGeometry[detID]) {
      return;
    }
    for (int iSt = 0; iSt < pIfs->GetNofTrackingStations(); ++iSt) {
      const auto& actVol  = pIfs->GetActiveVolume(iSt);
      const auto& fullVol = pIfs->GetFullVolume(iSt);
      fBuilder.AddLayer({
        .detId     = static_cast<int>(detID),              //  detector ID (ca::EDetectorID)
        .locId     = iSt,                                  //  local station ID
        .zRef      = actVol.GetCenterZ(),                  //  ref. z-position [cm]
        .zMin      = fullVol.GetMinZ(),                    //  lower bound in z-direction [cm]
        .zMax      = fullVol.GetMaxZ(),                    //  upper bound in z-direction [cm]
        .xMinFull  = fullVol.GetMinX(),                    //  lower bound of full volume in x-direction [cm]
        .xMaxFull  = fullVol.GetMaxX(),                    //  upper bound of full volume in x-direction [cm]
        .yMinFull  = fullVol.GetMinY(),                    //  lower bound of full volume in y-direction [cm]
        .yMaxFull  = fullVol.GetMaxY(),                    //  upper bound of full volume in y-direction [cm]
        .xMinAct   = actVol.GetMinX(),                     //  lower bound of active volume in x-direction [cm]
        .xMaxAct   = actVol.GetMaxX(),                     //  lower bound of active volume in x-direction [cm]
        .yMinAct   = actVol.GetMinY(),                     //  lower bound of active volume in x-direction [cm]
        .yMaxAct   = actVol.GetMaxY(),                     //  lower bound of active volume in x-direction [cm]
        .fieldType = GuessFieldType(actVol.GetCenterZ()),  //  field type
        .timeInfo  = pIfs->IsTimeInfoProvided(iSt)         //  time measurement flag
      });
    }
  };
 
  const auto& recoSetup = RecoSetupManager::Instance()->GetSetup();
  CollectStations(recoSetup.GetMvd(), EDetectorID::kMvd);
  CollectStations(recoSetup.GetSts(), EDetectorID::kSts);
  CollectStations(recoSetup.GetMuch(), EDetectorID::kMuch);
  CollectStations(recoSetup.GetTrd(), EDetectorID::kTrd);
  CollectStations(recoSetup.GetTof(), EDetectorID::kTof);
 
  // Set target properties
  const auto* pTarget = cbm::kf::Target::Instance();
  fBuilder.SetTargetProperty(pTarget->GetX(), pTarget->GetY(), pTarget->GetZ(), pTarget->GetDz(), pTarget->GetRmax());
 
  // Set material map creator
  auto pMaterialFactory = std::make_shared<MaterialMapFactory>();
  pMaterialFactory->SetSafeMaterialInitialization(kMatCreatorSafeMode);
  pMaterialFactory->SetDoRadialProjection(pTarget->GetZ());
  pMaterialFactory->SetNraysPerDim(kMatCreatorNrays);
  fBuilder.SetMaterialMapFactory(pMaterialFactory);
  fBuilder.SetMaterialCacheFile(MaterialCacheFile(), CbmSetup::Instance()->GetHash());
  fBuilder.Init();
  LOG(info) << "cbm::kf::TrackingGeoSetupFactory: initializing ... done";
}
catch (const std::exception& err) {
  LOG(fatal) << "cbm::kf::TrackingGeoSetupFactory: initializing ... failed. Reason: " << err.what();
}
 
 
// ---------------------------------------------------------------------------------------------------------------------
//
void TrackingGeoSetupFactory::CheckDetectorPresence(bool requireHitPresence)
{
  LOG(info) << "TrackingGeoSetupFactory: detector subsystems in geometry: ";
  auto Check = [&](EDetectorID detID) {
    auto modId              = cbm::algo::ca::ToCbmModuleId(detID);
    fvbDetInGeometry[detID] = CbmSetup::Instance()->IsActive(modId);
    if (requireHitPresence) {
      fvbDetHasHits[detID] = fvbDetInGeometry[detID] && FairRootManager::Instance()->GetObject(kDetHitBrName[detID]);
      LOG(info) << fmt::format("\t{:6} in geometry: {:5}, has hits: {:5}", std::string(cbm::util::ToString(modId)),
                               fvbDetInGeometry[detID], fvbDetHasHits[detID]);
    }
    else {
      fvbDetHasHits[detID] = fvbDetInGeometry[detID];
      LOG(info) << fmt::format("\t{:6} in geometry: {:5}", std::string(cbm::util::ToString(modId)),
                               fvbDetInGeometry[detID]);
    }
  };
 
  Check(EDetectorID::kMvd);
  Check(EDetectorID::kSts);
  Check(EDetectorID::kMuch);
  Check(EDetectorID::kTrd);
  Check(EDetectorID::kTof);
}
 
// ---------------------------------------------------------------------------------------------------------------------
//
std::string TrackingGeoSetupFactory::MaterialCacheFile() const
{
  namespace fs = boost::filesystem;
 
  auto setupTag = CbmSetup::Instance()->GetProvider()->GetSetup().GetTag();
  if (setupTag.empty()) {
    throw std::logic_error("The setup tag in CbmSetup is not defined");
  }
 
  TString sinkName     = FairRootManager::Instance()->GetSink()->GetFileName();
  auto sinkPath        = fs::path{sinkName.Data()};
  std::string sDataDir = sinkPath.parent_path().string();
  if (sDataDir.empty()) {
    sDataDir = ".";
  }
 
  return Form("%s/%s.mat.kf.bin", sDataDir.c_str(), setupTag.c_str());
}