diff --git a/algo/ca/core/pars/CaMaterialMonitor.h b/algo/ca/core/pars/CaMaterialMonitor.h index 82557f3286bf167920332a443bb96b955b9251ab..c220eef13db6f4e2299dcfb9953f90bceebabb65 100644 --- a/algo/ca/core/pars/CaMaterialMonitor.h +++ b/algo/ca/core/pars/CaMaterialMonitor.h @@ -83,7 +83,6 @@ namespace cbm::algo::ca /// get average radiation thickness among all passive bins double GetPassiveRadThickMean() const { return fPassiveRadThickMean; } - /// get the ration of hits that show up outside the material map double GetRatioOfOutsideHits() const { return fNhitsOutside / (fNhitsTotal + 1.e-8); } diff --git a/algo/kf/CMakeLists.txt b/algo/kf/CMakeLists.txt index ad6d4787c23e301370f65bac17e7904ad27d4a22..adc699bb9489a8009a46a43d05f011fec8b2d135 100644 --- a/algo/kf/CMakeLists.txt +++ b/algo/kf/CMakeLists.txt @@ -1 +1,2 @@ add_subdirectory(core) +add_subdirectory(tools) diff --git a/algo/kf/core/KfDefs.h b/algo/kf/core/KfDefs.h index 612bc9ebe0f99181f055eb551ac6ecd706b86b9c..abed6a7468806ba882536440240e55162b3d3c4b 100644 --- a/algo/kf/core/KfDefs.h +++ b/algo/kf/core/KfDefs.h @@ -31,13 +31,12 @@ namespace cbm::algo::kf /// \enum EFieldType /// \brief Magnetic field type in different setup regions - // TODO: SZh 21.08.2024: Develope the concept of field types at different setup regions + // TODO: SZh 21.08.2024: Develope the concept of field types at different setup regions. At the moment the only + // two global options are available. enum class EFieldType { - Target, ///< Field near the target - Tracker, ///< Field near the tracker subsystem - Low, ///< Magnetic field in the intermediate region (from tracker to null) - Null ///< No magnetic field + Normal, ///< Field near the tracker subsystem + Null ///< No magnetic field }; /// \struct Literal @@ -69,6 +68,7 @@ namespace cbm::algo::kf /// \brief Magnetic field function type /// \note The field function follows the FairField::GetFieldFunction signature using FieldFnFair_t = std::function<void(const double (&r)[3], double*)>; + } // namespace cbm::algo::kf namespace cbm::algo::kf::defs diff --git a/algo/kf/core/geo/KfSetup.h b/algo/kf/core/geo/KfSetup.h index fc26c6ac489bfa927439759f7ef4302a18a1aa37..7cced174893dfce18170b7ba611abaecd15b0d78 100644 --- a/algo/kf/core/geo/KfSetup.h +++ b/algo/kf/core/geo/KfSetup.h @@ -160,7 +160,7 @@ namespace cbm::algo::kf SetupInitializer(SetupInitializer&&) = delete; /// \brief Destructor - ~SetupInitializer() = default; + virtual ~SetupInitializer() = default; /// \brief Copy assignment operator SetupInitializer& operator=(const SetupInitializer&) = delete; @@ -181,6 +181,13 @@ namespace cbm::algo::kf /// \param material Material layer void AddMaterial(const MaterialMap& material); + /// \brief Virtual function, which initializes a KF-setup + /// \return true Initialization was successful + /// \return false Initialization failed + /// + /// This function can be re-implemented for a particular use-case... + virtual bool InitProperties() { return false; } + /// \brief Creates a setup instance /// \param bProvideOrigField Flag: true - defines optional variable of the original field template<typename T> @@ -227,6 +234,7 @@ namespace cbm::algo::kf std::set<MaterialMap> fMaterialLayers{}; ///< Set of material maps FieldFactory fFieldFactory; ///< Instance of field factory Target<double> fTarget; ///< Field to keep target properties + bool fbPropertiesInitialized{false}; ///< Flag: if all the necessary fields are initialized }; diff --git a/algo/kf/core/geo/KfTarget.h b/algo/kf/core/geo/KfTarget.h index 78008cc4501d81e6503ddf2c45d82dcccfb2c97e..60734745845fab50562d28b8841e832a7657aa01 100644 --- a/algo/kf/core/geo/KfTarget.h +++ b/algo/kf/core/geo/KfTarget.h @@ -112,9 +112,9 @@ namespace cbm::algo::kf ar& fZ; } - MaterialMap fMaterial; ///< Material map in the target region - T fX = defs::Undef<T>; ///< x-coordinate of the nominal target center [cm] - T fY = defs::Undef<T>; ///< y-coordinate of the nominal target center [cm] - T fZ = defs::Undef<T>; ///< z-coordinate of the nominal target center [cm] + MaterialMap fMaterial{}; ///< Material map in the target region + T fX{defs::Undef<T>}; ///< x-coordinate of the nominal target center [cm] + T fY{defs::Undef<T>}; ///< y-coordinate of the nominal target center [cm] + T fZ{defs::Undef<T>}; ///< z-coordinate of the nominal target center [cm] }; } // namespace cbm::algo::kf diff --git a/reco/kf/CMakeLists.txt b/reco/kf/CMakeLists.txt index 9a0ce4fe858b1bf06eb37a931d792b8e32d1d8a4..0bc036a5d297f3619c0ae68f70c68776414cf240 100644 --- a/reco/kf/CMakeLists.txt +++ b/reco/kf/CMakeLists.txt @@ -3,7 +3,7 @@ set(INCLUDE_DIRECTORIES ) set(SRCS - ${CMAKE_CURRENT_SOURCE_DIR}/CbmKfSetupInitializer.cxx + ${CMAKE_CURRENT_SOURCE_DIR}/CbmKfTrackingSetupInitializer.cxx ) SET_SOURCE_FILES_PROPERTIES(${SRCS} PROPERTIES COMPILE_FLAGS "-O3") @@ -41,7 +41,8 @@ target_link_libraries(CbmKf ROOT::Graf ROOT::Physics fmt::fmt - PRIVATE CbmRecoSts + PRIVATE KfTools + CbmRecoSts CbmSimSteer CbmRecoBase KFParticle @@ -67,7 +68,7 @@ install(TARGETS CbmKf DESTINATION lib) install( FILES - CbmKfSetupInitializer.h + CbmKfTrackingSetupInitializer.h DESTINATION include/ ) diff --git a/reco/kf/CbmKfSetupInitializer.cxx b/reco/kf/CbmKfSetupInitializer.cxx deleted file mode 100644 index ceac2c8c7d388a7bdc9d4e1bdf1b86182597c701..0000000000000000000000000000000000000000 --- a/reco/kf/CbmKfSetupInitializer.cxx +++ /dev/null @@ -1,16 +0,0 @@ -/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt - SPDX-License-Identifier: GPL-3.0-only - Authors: Sergei Zharko [committer] */ - -/// @file CbmKfSetupInitializer.h -/// @brief CBM-specific kf::Setup initializer (source) -/// @since 28.08.2024 -/// @author Sergei Zharko <s.zharko@gsi.de> - -#include "CbmKfSetupInitializer.h" - -using cbm::kf::SetupInitializer; - -// --------------------------------------------------------------------------------------------------------------------- -// -void SetupInitializer::Init() {} diff --git a/reco/kf/CbmKfSetupInitializer.h b/reco/kf/CbmKfSetupInitializer.h deleted file mode 100644 index 6b5489fbbf6bc0ef536e4f92bc7de44d33cfd6c8..0000000000000000000000000000000000000000 --- a/reco/kf/CbmKfSetupInitializer.h +++ /dev/null @@ -1,64 +0,0 @@ -/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt - SPDX-License-Identifier: GPL-3.0-only - Authors: Sergei Zharko [committer] */ - -/// @file CbmKfSetupInitializer.h -/// @brief CBM-specific kf::Setup initializer (header) -/// @since 28.08.2024 -/// @author Sergei Zharko <s.zharko@gsi.de> - -#pragma once - -#include "KfSetup.h" - -namespace cbm::kf -{ - /// \class SetupInitializer - /// \brief Encapsulation of the kf::Setup initialization routines for CBM - class SetupInitializer { - public: - /// \brief Default constructor - SetupInitializer() = default; - - /// \brief Copy constructor - SetupInitializer(const SetupInitializer&) = delete; - - /// \brief Move constructor - SetupInitializer(SetupInitializer&&) = delete; - - /// \brief Destructor - ~SetupInitializer() = default; - - /// \brief Copy assignment operator - SetupInitializer& operator=(const SetupInitializer&) = delete; - - /// \brief Move assignment operator - SetupInitializer& operator=(SetupInitializer&&) = delete; - - /// \brief Initializes the instance - /// \throws std::logic_error If initialization fails - void Init(); - - /// \brief Creates a kf::Setup instance - /// \tparam Underlying floating-point data-type for the created setup - template<typename T> - cbm::algo::kf::Setup<T> MakeSetup() const; - - private: - cbm::algo::kf::SetupInitializer fBaseInitializer; ///< Instance of the KF-core initializer - }; - - - // ********************************* - // ** Template methods definition ** - // ********************************* - - // ------------------------------------------------------------------------------------------------------------------- - // - template<typename T> - cbm::algo::kf::Setup<T> cbm::kf::SetupInitializer::MakeSetup() const - { - using cbm::algo::kf::Setup; - return Setup<T>(); - } -} // namespace cbm::kf diff --git a/reco/kf/CbmKfTrackingSetupInitializer.cxx b/reco/kf/CbmKfTrackingSetupInitializer.cxx new file mode 100644 index 0000000000000000000000000000000000000000..bb6f465d3a1f0a1307093f487822b5d0055259a2 --- /dev/null +++ b/reco/kf/CbmKfTrackingSetupInitializer.cxx @@ -0,0 +1,200 @@ +/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt + SPDX-License-Identifier: GPL-3.0-only + Authors: Sergei Zharko [committer] */ + +/// \file CbmKfTrackingSetupInitializer.cxx +/// \brief Tracking setup initializer in CBM (source) +/// \since 28.08.2024 +/// \author Sergei Zharko <s.zharko@gsi.de> + +#include "CbmKfTrackingSetupInitializer.h" + +#include "CbmMuchTrackingInterface.h" +#include "CbmMvdTrackingInterface.h" +#include "CbmSetup.h" +#include "CbmStsTrackingInterface.h" +#include "CbmTofTrackingInterface.h" +#include "CbmTrdTrackingInterface.h" +#include "FairField.h" +#include "FairRunAna.h" +#include "KfDefs.h" +#include "KfMaterialMapCreator.h" +#include "TGeoManager.h" +#include "TGeoNode.h" +#include "TGeoTube.h" +#include "TGeoVolume.h" + +#include <functional> + +using cbm::kf::TrackingSetupInitializer; +using kf::tools::MaterialMapCreator; + +// --------------------------------------------------------------------------------------------------------------------- +// +bool TrackingSetupInitializer::InitProperties() +try { + using cbm::algo::kf::EFieldType; + + // Check, if a subsystem exists in the setup + fbUseMvd &= CbmSetup::Instance()->IsActive(ECbmModuleId::kMvd); + fbUseSts &= CbmSetup::Instance()->IsActive(ECbmModuleId::kSts); + fbUseMuch &= CbmSetup::Instance()->IsActive(ECbmModuleId::kMuch); + fbUseTrd &= CbmSetup::Instance()->IsActive(ECbmModuleId::kTrd); + fbUseTof &= CbmSetup::Instance()->IsActive(ECbmModuleId::kTof); + + // Magnetic field initialization + if (FairRunAna::Instance()->GetField()) { + this->SetFieldFunction( + [](const double(&r)[3], double* B) { FairRunAna::Instance()->GetField()->GetFieldValue(r, B); }, + EFieldType::Normal); + } + else { + this->SetFieldFunction(cbm::algo::kf::defs::ZeroFieldFn, EFieldType::Null); + } + + // Tracking station property initialization + { + // Geometrical layer (contains information on the station sizes) + // + struct GeoLayer { + double zRef; // Reference z-coordinate [cm] + double zMin; // Lower z-coordinate boundary [cm] + double zMax; // Upper z-coordinate boundary [cm] + double xMax; // Half size of the layer in the x-direction [cm] + double yMax; // Half size of the layer in the y-direction [cm] + bool operator<(const GeoLayer& rhs) const { return zRef < rhs.zRef; } + }; + std::set<GeoLayer> geoStations; + + auto CollectStations = [&](const auto* pIfs) -> void { + if (!pIfs) { + return; + } + for (int iSt = 0; iSt < pIfs->GetNtrackingStations(); ++iSt) { + geoStations.emplace(GeoLayer{.zRef = pIfs->GetZref(iSt), + .zMin = pIfs->GetZmin(iSt), + .zMax = pIfs->GetZmax(iSt), + .xMax = std::max(std::fabs(pIfs->GetXmin(iSt)), std::fabs(pIfs->GetXmax(iSt))), + .yMax = std::max(std::fabs(pIfs->GetYmin(iSt)), std::fabs(pIfs->GetYmax(iSt)))}); + } + }; + CollectStations(fbUseMvd ? CbmMvdTrackingInterface::Instance() : nullptr); + CollectStations(fbUseSts ? CbmStsTrackingInterface::Instance() : nullptr); + CollectStations(fbUseMuch ? CbmMuchTrackingInterface::Instance() : nullptr); + CollectStations(fbUseTrd ? CbmTrdTrackingInterface::Instance() : nullptr); + CollectStations(fbUseTof ? CbmTofTrackingInterface::Instance() : nullptr); + + // Retriev target properties + // TODO: Provide some beter handler for target properties access + auto [targX, targY, targZ, targHalfThick, targRadius] = this->GetTarget(); + + // Init material map creator + MaterialMapCreator materialCreator{}; + materialCreator.SetSafeMaterialInitialization(kMatCreatorSafeMode); + materialCreator.SetDoRadialProjection(targZ); + materialCreator.SetNraysPerDim(kMatCreatorNrays); + + // Target initialization + double zLast = targZ + targHalfThick + kTargMaterialOffset; // z-coordinate of the target material upper limit + { + // Material of the target + double xyMax{1.3 * targRadius}; + int nBins{std::min(static_cast<int>(std::ceil(2. * xyMax / kMatCreatorPitch)), kMatCreatorMaxNbins)}; + assert(nBins > 0); + double zFirst{targZ - targHalfThick}; // z-coordinate of the target material lower limit + LOG(info) << "- collecting target material (z = " << zFirst << " - " << zLast << ", size = " << xyMax << ")\n"; + this->SetTargetProperties(targX, targY, targZ, kTargFieldInitStep, + materialCreator.GenerateMaterialMap(targZ, zFirst, zLast, xyMax, nBins)); + } + + // Material and field layers initialization + for (auto layerIt = geoStations.begin(); layerIt != geoStations.end(); ++layerIt) { + double z1 = layerIt->zMax; + double z2 = z1; + if (std::next(layerIt) != geoStations.end()) { + z2 = std::next(layerIt)->zMin; + } + double zNew = 0.5 * (z1 + z2); + double xyMax = 1.3 * std::max(std::fabs(layerIt->xMax), std::fabs(layerIt->yMax)); + int nBins = std::min(static_cast<int>(std::ceil(2. * xyMax / kMatCreatorPitch)), kMatCreatorMaxNbins); + assert(nBins > 0); + + LOG(info) << "- collecting material between z = " << zLast << " and " << zNew << " [cm]"; + this->AddMaterial(materialCreator.GenerateMaterialMap(layerIt->zRef, zLast, zNew, xyMax, nBins)); + this->AddFieldSliceReference(layerIt->xMax, layerIt->yMax, layerIt->zRef); + zLast = zNew; + } + } + return true; +} +catch (const std::exception& err) { + LOG(error) << "kf::TrackingSetupInitializer: Tracking setup was not initialized. Reason: " << err.what(); + return false; +} + +// --------------------------------------------------------------------------------------------------------------------- +// +void TrackingSetupInitializer::Use(bool mvd, bool sts, bool much, bool trd, bool tof) +{ + fbUseMvd = mvd; + fbUseSts = sts; + fbUseMuch = much; + fbUseTrd = trd; + fbUseTof = tof; +} + +// --------------------------------------------------------------------------------------------------------------------- +// +std::tuple<double, double, double, double, double> TrackingSetupInitializer::GetTarget() +{ + std::function<void(TString&, TGeoNode*)> FindTargetNode = [&](TString& targetPath, TGeoNode* targetNode) -> void { + if (!targetNode) { // init at the top of the tree + targetNode = gGeoManager->GetTopNode(); + targetPath = "/" + TString(targetNode->GetName()); + } + + if (TString(targetNode->GetName()).Contains("target")) { + return; + } + + for (Int_t iNode = 0; iNode < targetNode->GetNdaughters(); ++iNode) { + TGeoNode* newNode = targetNode->GetDaughter(iNode); + TString newPath = targetPath + "/" + newNode->GetName(); + FindTargetNode(newPath, newNode); + if (newNode) { + targetPath = newPath; + targetNode = newNode; + return; + } + } + targetPath = ""; + targetNode = nullptr; + }; + + TString targetPath; + TGeoNode* pTargetNode{nullptr}; + FindTargetNode(targetPath, pTargetNode); + + if (!pTargetNode) { + throw std::runtime_error("kf::TrackingSetupInitializer: target node is not found in the setup"); + } + + Double_t local[3] = {0., 0., 0.}; // target centre, local c.s. + Double_t global[3]; // target centre, global c.s. + gGeoManager->cd(targetPath); + gGeoManager->GetCurrentMatrix()->LocalToMaster(local, global); + + double halfThick{cbm::algo::kf::defs::Undef<double>}; + double outerRadius{cbm::algo::kf::defs::Undef<double>}; + if (const auto* pTube = dynamic_cast<TGeoTube*>(pTargetNode->GetVolume()->GetShape())) { + halfThick = pTube->GetDz(); + outerRadius = pTube->GetRmax(); + } + else { + throw std::logic_error("kf::TrackingSetupInitializer: target is supposed to be a Tube, but it is not. Please, " + "provide a proper handling of the new target shape (return it's reference central point " + "and half of its thickness)"); + } + + return std::make_tuple(local[0], local[1], local[2], halfThick, outerRadius); +} diff --git a/reco/kf/CbmKfTrackingSetupInitializer.h b/reco/kf/CbmKfTrackingSetupInitializer.h new file mode 100644 index 0000000000000000000000000000000000000000..62f22cb68dd86bdd651eefee671ea0ea961f97ed --- /dev/null +++ b/reco/kf/CbmKfTrackingSetupInitializer.h @@ -0,0 +1,56 @@ +/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt + SPDX-License-Identifier: GPL-3.0-only + Authors: Sergei Zharko [committer] */ + +/// \file CbmKfTrackingSetupInitializer.h +/// \brief Tracking setup initializer in CBM (source) +/// \since 28.08.2024 +/// \author Sergei Zharko <s.zharko@gsi.de> + +#pragma once + +#include "KfSetup.h" + +#include <tuple> + +namespace cbm::kf +{ + /// \class TrackingSetupInitializer + /// \brief Encapsulation of the kf::Setup initialization routines for CBM + class TrackingSetupInitializer : public cbm::algo::kf::SetupInitializer { + public: + /// \brief Default constructor + TrackingSetupInitializer(); + + /// \brief Initializes the instance + bool InitProperties() override; + + /// \brief Enables/disables detector subsystems in the setup + /// \param mvd Is MVD used + /// \param sts Is STS used + /// \param much Is MuCh used + /// \param trd Is TRD used + /// \param tof Is TOF used + void Use(bool mvd, bool sts, bool much, bool trd, bool tof); + + private: + /// \brief Gets target position and half-thickness + /// \return tuple(x, y, z, half-thickness, outer radius) [cm] + std::tuple<double, double, double, double, double> GetTarget(); + + // Material map creator properties (TODO: Provide setters, if needed) + static constexpr double kMatCreatorPitch{0.1}; ///< Material budget map minimal bin size [cm] + static constexpr int kMatCreatorMaxNbins{100}; ///< Max number of bins in the material budget map in x(y) axis + static constexpr int kMatCreatorNrays{3}; ///< Number of rays per dimension for the material budget + static constexpr bool kMatCreatorSafeMode{true}; ///< Safe mode of the material map creation + static constexpr double kTargFieldInitStep{2.5}; ///< Step between nodes in the target field initialization [cm] + static constexpr double kTargMaterialOffset{ + 0.2}; ///< Offset between target volume and the corresponding material zMax [cm] + + bool fbUseMvd{false}; ///< Are MVD stations included in the tracking setup + bool fbUseSts{false}; ///< Are STS stations included in the tracking setup + bool fbUseMuch{false}; ///< Are MuCh stations included in the tracking setup + bool fbUseTrd{false}; ///< Are TRD stations included in the tracking setup + bool fbUseTof{false}; ///< Are TOF stations included in the tracking setup + }; +} // namespace cbm::kf