diff --git a/algo/ca/core/CMakeLists.txt b/algo/ca/core/CMakeLists.txt
index e54aeda0763ca94e149c3cd664e892987ceec700..2be924997474f03ae9311ec060bd857d21cf725a 100644
--- a/algo/ca/core/CMakeLists.txt
+++ b/algo/ca/core/CMakeLists.txt
@@ -19,7 +19,6 @@ set(SRCS
${CMAKE_CURRENT_SOURCE_DIR}/pars/CaConfigReader.cxx
${CMAKE_CURRENT_SOURCE_DIR}/pars/CaInitManager.cxx
${CMAKE_CURRENT_SOURCE_DIR}/pars/CaIteration.cxx
- ${CMAKE_CURRENT_SOURCE_DIR}/pars/CaMaterialMap.cxx
${CMAKE_CURRENT_SOURCE_DIR}/pars/CaParameters.cxx
${CMAKE_CURRENT_SOURCE_DIR}/pars/CaSearchWindow.cxx
${CMAKE_CURRENT_SOURCE_DIR}/pars/CaStation.cxx
@@ -112,7 +111,6 @@ install(
pars/CaDefs.h
pars/CaInitManager.h
pars/CaIteration.h
- pars/CaMaterialMap.h
pars/CaParameters.h
pars/CaSearchWindow.h
pars/CaStation.h
diff --git a/algo/ca/core/pars/CaInitManager.cxx b/algo/ca/core/pars/CaInitManager.cxx
index 8edfa2a2c5f9e37a6fa217ba0ead02419fcad6ac..02530363ecbd5dff79a54d9ba886182f6cd857bb 100644
--- a/algo/ca/core/pars/CaInitManager.cxx
+++ b/algo/ca/core/pars/CaInitManager.cxx
@@ -23,7 +23,6 @@
using namespace constants;
using cbm::algo::ca::InitManager;
-using cbm::algo::ca::MaterialMap;
using cbm::algo::ca::Station;
using cbm::algo::ca::StationInitializer;
diff --git a/algo/ca/core/pars/CaMaterialMap.cxx b/algo/ca/core/pars/CaMaterialMap.cxx
deleted file mode 100644
index 75fc25c5721476688d9c3e0ec80f94db867def79..0000000000000000000000000000000000000000
--- a/algo/ca/core/pars/CaMaterialMap.cxx
+++ /dev/null
@@ -1,146 +0,0 @@
-/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
- SPDX-License-Identifier: GPL-3.0-only
- Authors: Sergey Gorbunov, Sergei Zharko [committer] */
-
-#include "CaMaterialMap.h"
-
-#include "AlgoFairloggerCompat.h"
-
-#include <iomanip>
-#include <sstream>
-#include <vector>
-
-/********************
- * MaterialMap class *
- ********************/
-
-using namespace cbm::algo::ca;
-
-//------------------------------------------------------------------------------------------------------------------------------------
-//
-MaterialMap::MaterialMap(MaterialMap&& other) noexcept { this->Swap(other); }
-
-//------------------------------------------------------------------------------------------------------------------------------------
-//
-MaterialMap& MaterialMap::operator=(MaterialMap&& other) noexcept
-{
- if (this != &other) {
- MaterialMap tmp(std::move(other));
- this->Swap(tmp);
- }
- return *this;
-}
-
-//------------------------------------------------------------------------------------------------------------------------------------
-int MaterialMap::GetBin(float x, float y) const
-{
- int i = static_cast<int>((x + fXYmax) * fFactor);
- int j = static_cast<int>((y + fXYmax) * fFactor);
- if (i < 0 || j < 0 || i >= fNbins || j >= fNbins) {
- return -1;
- }
- return i + j * fNbins;
-}
-
-//------------------------------------------------------------------------------------------------------------------------------------
-//
-float MaterialMap::GetRadThickScal(float x, float y) const
-{
- x = (x < fXYmax && x >= -fXYmax) ? x : 0;
- y = (y < fXYmax && y >= -fXYmax) ? y : 0;
- int i = static_cast<int>((x + fXYmax) * fFactor);
- int j = static_cast<int>((y + fXYmax) * fFactor);
- i = (i < fNbins && i >= 0) ? i : fNbins / 2;
- j = (j < fNbins && j >= 0) ? j : fNbins / 2;
- return fTable[i + j * fNbins];
-}
-
-//------------------------------------------------------------------------------------------------------------------------------------
-//
-fvec MaterialMap::GetRadThickVec(fvec x, fvec y) const
-{
- fvec r;
- for (size_t i = 0; i < fvec::size(); i++)
- r[i] = GetRadThickScal(x[i], y[i]);
- return r;
-}
-
-//------------------------------------------------------------------------------------------------------------------------------------
-//
-void MaterialMap::CheckConsistency() const
-{
- /* (i) Check, if the object was initialized */
- if (IsNaN()) {
- throw std::logic_error("MaterialMap: class was not initialized");
- }
-
- /* (ii) Check if the thickness values correct (non-negative) */
- bool isThicknessOk = true;
- for (int i = 0; i < int(fTable.size()); ++i) {
- if (fTable[i] < 0) {
- isThicknessOk = false;
- LOG(error) << "MaterialMap: found illegal negative station thickness value " << fTable[i]
- << " at iBinX = " << (i % fNbins) << ", iBin = " << (i / fNbins);
- }
- }
- if (!isThicknessOk) {
- throw std::logic_error("MaterialMap: incorrect station thickness values found in the thickness map");
- }
-}
-
-//------------------------------------------------------------------------------------------------------------------------------------
-//
-void MaterialMap::Initialize(int nBins, float xyMax, float zRef, float zMin, float zMax)
-{
- fNbins = nBins;
- fXYmax = xyMax;
- fZref = zRef;
- fZmin = zMin;
- fZmax = zMax;
-
- if (fNbins < 1) {
- std::stringstream aStream;
- aStream << "MaterialMap: object cannot be initialized with non-positive nBins = " << fNbins;
- throw std::logic_error(aStream.str());
- }
-
- if (fXYmax < 0.) {
- std::stringstream aStream;
- aStream << "MaterialMap: object cannot be initialized with non-positive XYmax = " << fXYmax << " [cm]";
- throw std::logic_error(aStream.str());
- }
-
- if (!((fZmin <= fZref) && (fZref <= zMax))) {
- std::stringstream aStream;
- aStream << "MaterialMap: object cannot be initialized with inconsistent Z: min " << fZmin << " ref " << fZref
- << " max " << fZmax << " [cm]";
- throw std::logic_error(aStream.str());
- }
-
- fFactor = 0.5 * fNbins / fXYmax;
- fTable.resize(fNbins * fNbins);
-}
-
-//------------------------------------------------------------------------------------------------------------------------------------
-//
-void MaterialMap::Swap(MaterialMap& other) noexcept
-{
- std::swap(fNbins, other.fNbins);
- std::swap(fXYmax, other.fXYmax);
- std::swap(fFactor, other.fFactor);
- std::swap(fZref, other.fZref);
- std::swap(fZmin, other.fZmin);
- std::swap(fZmax, other.fZmax);
- std::swap(fTable, other.fTable); // Probably can cause segmentation violation (did not understand)
-}
-
-//------------------------------------------------------------------------------------------------------------------------------------
-//
-std::string MaterialMap::ToString() const
-{
- std::stringstream msg;
- using std::setw;
- msg << "[mat.map] nBins " << setw(6) << fNbins << ", XYmax " << setw(12) << fXYmax << " , z (ref, min, max) "
- << setw(12) << fZref << ' ' << setw(12) << fZmin << ' ' << setw(12) << fZmax;
- return msg.str();
-}
diff --git a/algo/ca/core/pars/CaMaterialMap.h b/algo/ca/core/pars/CaMaterialMap.h
deleted file mode 100644
index f72ebb8d31dacda8ffec447dd86851b69a37a938..0000000000000000000000000000000000000000
--- a/algo/ca/core/pars/CaMaterialMap.h
+++ /dev/null
@@ -1,136 +0,0 @@
-/* Copyright (C) 2007-2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
- SPDX-License-Identifier: GPL-3.0-only
- Authors: Igor Kulakov, Sergey Gorbunov [committer], Andrey Lebedev, Sergei Zharko */
-
-#pragma once // include this header only once per compilation unit
-
-#include "CaSimd.h"
-
-#include <boost/serialization/vector.hpp>
-
-#include <iomanip>
-#include <sstream>
-#include <string>
-#include <vector>
-
-namespace cbm::algo::ca
-{
- /// \class MaterialMap
- /// \brief A map of station thickness in units of radiation length (X0) to the specific point in XY plane
- class MaterialMap {
- public:
- /// \brief Default constructor
- MaterialMap() = default;
-
- /// \brief Copy constructor
- MaterialMap(const MaterialMap& other) = default;
-
- /// \brief Copy assignment operator
- MaterialMap& operator=(const MaterialMap& other) = default;
-
- /// \brief Move constructor
- MaterialMap(MaterialMap&& other) noexcept;
-
- /// \brief Move assignment operator
- MaterialMap& operator=(MaterialMap&& other) noexcept;
-
- /// \brief Destructor
- ~MaterialMap() noexcept = default;
-
- /// \brief Gets number of bins (rows or columns) of the material table
- int GetNbins() const { return fNbins; }
-
- /// \brief Gets radius in cm of the material table
- float GetXYmax() const { return fXYmax; }
-
- /// \brief Gets reference Z of the material in cm
- float GetZref() const { return fZref; }
-
- /// \brief Gets minimal Z of the collected material in cm
- float GetZmin() const { return fZmin; }
-
- /// \brief Gets maximal Z of the collected material in cm
- float GetZmax() const { return fZmax; }
-
- /// \brief Gets value of X/X0 in a given cell of the material table by the indeces of the row and column
- /// \param iBinX Index of table column
- /// \param iBinY Index of table row
- float GetRadThickBin(int iBinX, int iBinY) const { return fTable[iBinX + fNbins * iBinY]; }
-
- /// \brief Gets value of X/X0 in a given cell of the material table by the index of the bin
- /// \param iBin Index of the bin in 2d table
- float GetRadThickBin(int iBin) const { return fTable[iBin]; }
-
- /// \brief Gets material thickness in units of X0 in (x,y) point of the station
- /// \param x X coordinate of the point [cm]
- /// \param y Y coordinate of the point [cm]
- float GetRadThickScal(float x, float y) const;
-
- /// \brief Gets material thickness in units of X0 in (x,y) point of the station
- /// cbm::algo::ca::fvec type can be float, that is why "Vec" and "Scal" specifications
- /// \param x X coordinate of the point [cm] (SIMDized vector)
- /// \param y Y coordinate of the point [cm] (SIMDized veotor)
- fvec GetRadThickVec(fvec x, fvec y) const;
-
- /// \brief Checks, if the fields are NaN
- bool IsNaN() const
- {
- return kfutils::IsUndefined(fNbins) || kfutils::IsUndefined(fXYmax) || kfutils::IsUndefined(fFactor)
- || kfutils::IsUndefined(fZref) || kfutils::IsUndefined(fZmin) || kfutils::IsUndefined(fZmax);
- }
-
- /// \brief Verifies class invariant consistency
- void CheckConsistency() const;
-
- /// \brief Sets value of material thickness in units of X0 for a given cell of the material table
- /// \param iBinX Index of table column
- /// \param iBinY Index of table row
- /// \param thickness Thickness of the material in units of X0
- /// \note Indices of rows and columns in the table runs from 0 to nBins-1 inclusively, where nBins is the number
- /// both of rows and columns. One should be careful while reading and storing the table from ROOT-file,
- /// because iBinX = 0 and iBinY = 0 in the TH1::SetBinContent method of usually defines the underflow bin.
- void SetRadThickBin(int iBinX, int iBinY, float thickness) { fTable[iBinX + fNbins * iBinY] = thickness; }
-
- /// \brief Sets properties of the material table -- number of rows or columnts and the size of station in XY plane
- /// \param nBins Number of rows or columns
- /// \param xyMax Size of station in x and y dimensions [cm]
- void Initialize(int nBins, float xyMax, float zRef, float zMin, float zMax);
-
- /// \brief Swap method
- void Swap(MaterialMap& other) noexcept;
-
- /// \brief Get bin index for (x,y). Returns -1 when outside of the map
- int GetBin(float x, float y) const;
-
- /// \brief String representation of the object
- std::string ToString() const;
-
- /// \brief Comparison operator
- /// \note To be used for the material map ordering along the z-axis
- bool operator<(const MaterialMap& r) const { return fZref < r.fZref; }
-
- private:
- int fNbins = kfdefs::Undef<int>; ///< Number of rows (== N columns) in the material budget table
- float fXYmax = kfdefs::Undef<float>; ///< Size of the station in x and y dimensions [cm]
- float fFactor = kfdefs::Undef<float>; ///< Util. var. for the conversion of point coordinates to row/column id
- float fZref = kfdefs::Undef<float>; ///< Reference Z of the collected material [cm]
- float fZmin = kfdefs::Undef<float>; ///< Minimal Z of the collected material [cm]
- float fZmax = kfdefs::Undef<float>; ///< Minimal Z of the collected material [cm]
- std::vector<float> fTable{}; ///< Material budget table
-
- /// \brief Serialization function
- friend class boost::serialization::access;
- template<class Archive>
- void serialize(Archive& ar, const unsigned int)
- {
- ar& fNbins;
- ar& fXYmax;
- ar& fFactor;
- ar& fZref;
- ar& fZmin;
- ar& fZmax;
- ar& fTable;
- }
- } _fvecalignment;
-
-} // namespace cbm::algo::ca
diff --git a/algo/ca/core/pars/CaParameters.cxx b/algo/ca/core/pars/CaParameters.cxx
index 3561da863e5e5f78adfccd006e18d3b071bb36bb..c89279a12c1c24c77b9b03cc28d931f4ad47b931 100644
--- a/algo/ca/core/pars/CaParameters.cxx
+++ b/algo/ca/core/pars/CaParameters.cxx
@@ -15,7 +15,6 @@
using cbm::algo::ca::EDetectorID;
using cbm::algo::ca::Iteration;
-using cbm::algo::ca::MaterialMap;
using cbm::algo::ca::Parameters;
using cbm::algo::ca::kfutils::CheckSimdVectorEquality;
diff --git a/algo/ca/core/pars/CaParameters.h b/algo/ca/core/pars/CaParameters.h
index e2f145c0ed4affe124b7c8b4f940d18808be3c3d..dbbf456e6d73777445aeb9f6b7076a8fd3fef464 100644
--- a/algo/ca/core/pars/CaParameters.h
+++ b/algo/ca/core/pars/CaParameters.h
@@ -11,7 +11,6 @@
#include "CaDefs.h"
#include "CaIteration.h"
-#include "CaMaterialMap.h"
#include "CaSearchWindow.h"
#include "CaStation.h"
#include "CaVector.h"
@@ -34,10 +33,9 @@ namespace cbm::algo::ca
using constants::Undef;
/// Type definitions for used containers
- using IterationsContainer_t = cbm::algo::ca::Vector<cbm::algo::ca::Iteration>;
+ using IterationsContainer_t = Vector<Iteration>;
template<typename DataT>
using StationsContainer_t = std::array<ca::Station<DataT>, constants::size::MaxNstations>;
- using MaterialContainer_t = std::array<ca::MaterialMap, constants::size::MaxNstations>;
enum TrackingMode
{
diff --git a/algo/ca/core/pars/CaStationInitializer.cxx b/algo/ca/core/pars/CaStationInitializer.cxx
index 6c48c65c1781a26c171c7d53294f820043a43f8b..f99927d4a0dcba0d032fcc7082b916696568f939 100644
--- a/algo/ca/core/pars/CaStationInitializer.cxx
+++ b/algo/ca/core/pars/CaStationInitializer.cxx
@@ -19,7 +19,6 @@
using cbm::algo::ca::EDetectorID;
using cbm::algo::ca::fvec;
-using cbm::algo::ca::MaterialMap;
using cbm::algo::ca::Station;
using cbm::algo::ca::StationInitializer;
diff --git a/algo/ca/core/pars/CaStationInitializer.h b/algo/ca/core/pars/CaStationInitializer.h
index 7347cee1300ddcb2df15374312248ed1005f3ad1..3ac7fb349d7aa87108567c1a3faed3414db35219 100644
--- a/algo/ca/core/pars/CaStationInitializer.h
+++ b/algo/ca/core/pars/CaStationInitializer.h
@@ -9,7 +9,6 @@
#pragma once // include this header only once per compilation unit
-#include "CaMaterialMap.h"
#include "CaObjectInitController.h"
#include "CaSimd.h"
#include "CaStation.h"
diff --git a/reco/L1/CMakeLists.txt b/reco/L1/CMakeLists.txt
index 8c61bf88f11d7fdec12c422ab765c00b9b659dff..6f6b1a79af59ad98b7b49317f881c966106838f6 100644
--- a/reco/L1/CMakeLists.txt
+++ b/reco/L1/CMakeLists.txt
@@ -48,7 +48,6 @@ set(SRCS
catools/CaToolsDebugger.cxx
catools/CaToolsWindowFinder.cxx
catools/CaToolsWFExpression.cxx
- catools/CaToolsMaterialHelper.cxx
qa/CbmCaInputQaBase.cxx
qa/CbmCaInputQaMvd.cxx
@@ -83,7 +82,6 @@ set(HEADERS
catools/CaToolsHitRecord.h
catools/CaToolsDef.h
utils/CbmCaIdealHitProducer.h
- catools/CaToolsMaterialHelper.h
catools/CaToolsField.h
qa/CbmCaInputQaBase.h
qa/CbmCaHitQaData.h
diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx
index fd3eff0f43d8882418717c796460ee11136c7a50..a66957edd371d1ba14d94fcc58e34290e9a20be4 100644
--- a/reco/L1/CbmL1.cxx
+++ b/reco/L1/CbmL1.cxx
@@ -20,7 +20,6 @@
#include "CbmL1.h"
-#include "CaToolsMaterialHelper.h"
#include "CbmMCDataManager.h"
#include "CbmMuchTrackingInterface.h"
#include "CbmMvdTrackingInterface.h"
@@ -74,7 +73,6 @@ using CaTrack = cbm::algo::ca::Track;
using cbm::algo::ca::Parameters;
using cbm::ca::MCModule;
using cbm::ca::TimeSliceReader;
-using cbm::ca::tools::MaterialHelper;
ClassImp(CbmL1);
diff --git a/reco/L1/catools/CaToolsMaterialHelper.cxx b/reco/L1/catools/CaToolsMaterialHelper.cxx
deleted file mode 100644
index a12c1e450a3fdccf593848c3658cf14b74e40c05..0000000000000000000000000000000000000000
--- a/reco/L1/catools/CaToolsMaterialHelper.cxx
+++ /dev/null
@@ -1,404 +0,0 @@
-/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
- SPDX-License-Identifier: GPL-3.0-only
- Authors: Sergey Gorbunov [committer] */
-
-
-#include "CaToolsMaterialHelper.h"
-
-#include "TGeoManager.h"
-#include "TGeoMedium.h"
-#include "TGeoNavigator.h"
-#include "TGeoNode.h"
-#include "TGeoPhysicalNode.h"
-#include "TGeoVoxelFinder.h"
-
-#include <Logger.h>
-
-#include <sstream>
-#include <thread>
-
-using cbm::ca::tools::MaterialHelper;
-
-using cbm::algo::ca::MaterialMap;
-
-ClassImp(MaterialHelper);
-
-MaterialHelper::MaterialHelper()
-{
- // constructor
-
- // save the current number of threads to restore it when the work is done
- fNthreadsOld = gGeoManager->GetMaxThreads();
-
- // set defaut n threads to 1
-
- fNthreads = 1;
-
- LOG(info) << "------------------------------------------------------------------------------------";
- LOG(info) << " create L1 material maps ";
- LOG(info) << "";
- LOG(info) << " !! To run it with the full speed, set: \"export OMP_NUM_THREADS=<n CPUs>\" before running !!";
- LOG(info) << "";
- // if possible, read the allowed n threads from the environment variable
- {
- const char* env = std::getenv("OMP_NUM_THREADS");
- if (env) {
- fNthreads = TString(env).Atoi();
- LOG(info) << " found environment variable OMP_NUM_THREADS = \"" << env << "\", read as integer: " << fNthreads;
- }
- }
-
- // ensure that at least one thread is set
-
- if (fNthreads < 1) {
- fNthreads = 1;
- }
-
- LOG(info) << " Maps will be created using " << fNthreads << " CPU threads";
- LOG(info) << "";
- LOG(info) << " It might crash because of a ROOT bug. In this case do one of the following: ";
- LOG(info) << "";
- LOG(info) << " - fix illegal overlaps in the geometry via gGeoManager->CheckOverlaps()";
- LOG(info) << " - run with CbmL1::SetSafeMaterialInitialization() option. ";
- LOG(info) << " It fixes the crash, but slows down the initialization significantly.";
- LOG(info) << " - manually disable voxel finders for problematic volumes in CaToolsMaterialHelper.cxx";
- LOG(info) << "";
- LOG(info) << "------------------------------------------------------------------------------------";
-
- InitThreads();
-
- //
- // It looks like TGeoVoxelFinder contains a bug and therefore sometimes produces a crash.
- //
- // Currently, the crash only happens on some of TGeoVolumeAssembly volumes
- // and only when an (inproper?) alignment matrices are applied.
- //
- // Here we try to catch this case and disable those finders.
- //
- // Disabling all the voxel finders leads to a very long initialization, try to avoid it!
- //
- // TODO: make a bug report to the ROOT team; remove the following code once the bug is fixed.
- //
-
- if (fDoSafeInitialization) { // disable all voxel finders
- TObjArray* volumes = gGeoManager->GetListOfVolumes();
- if (volumes) {
- for (int iv = 0; iv < volumes->GetEntriesFast(); iv++) {
- TGeoVolume* vol = dynamic_cast<TGeoVolume*>(volumes->At(iv));
- if (!vol) {
- continue;
- }
- TGeoVoxelFinder* finder = vol->GetVoxels();
- if (finder) {
- finder->SetInvalid();
- }
- }
- }
- }
- else { // disable some voxel finders in some cases
-
- // check if any alignment was applied
-
- bool isAlignmentApplied = false;
- {
- TObjArray* nodes = gGeoManager->GetListOfPhysicalNodes();
- if (nodes) {
- for (int in = 0; in < nodes->GetEntriesFast(); in++) {
- TGeoPhysicalNode* node = dynamic_cast<TGeoPhysicalNode*>(nodes->At(in));
- if (!node) continue;
- if (node->IsAligned()) {
- isAlignmentApplied = true;
- break;
- }
- }
- }
- }
-
- // disable potentially problematic voxel finders
-
- if (isAlignmentApplied) {
- TObjArray* volumes = gGeoManager->GetListOfVolumes();
- if (volumes) {
- for (int iv = 0; iv < volumes->GetEntriesFast(); iv++) {
- TGeoVolumeAssembly* vol = dynamic_cast<TGeoVolumeAssembly*>(volumes->At(iv));
- if (!vol) {
- continue;
- }
- TGeoVoxelFinder* finder = vol->GetVoxels();
- if (finder) {
- finder->SetInvalid();
- }
- }
- }
- }
- }
-}
-
-
-MaterialHelper::~MaterialHelper()
-{
- // destructor
-
- CleanUpThreads();
-
- // delete the navigators, created in this class
-
- for (auto i = fNavigators.begin(); i != fNavigators.end(); i++) {
- gGeoManager->RemoveNavigator(*i);
- }
-
- // once TGeoManager is swithched in multithreaded mode, there is no way to make it non-mltithreaded again
- // therefore we should set SetMaxThreads( >=0 )
-
- if (fNthreadsOld > 0) {
- gGeoManager->SetMaxThreads(fNthreadsOld);
- }
- else {
- gGeoManager->SetMaxThreads(1);
- }
-
- // ensure that the default navigator exists
-
- GetCurrentNavigator(-1);
-}
-
-
-TGeoNavigator* MaterialHelper::GetCurrentNavigator(int iThread)
-{
- // Get navigator for current thread, create it if it doesnt exist.
- // Produce an error when anything goes wrong
-
- TGeoNavigator* navi = gGeoManager->GetCurrentNavigator();
- if (!navi) {
-
- navi = gGeoManager->AddNavigator();
-
- if (iThread >= 0) {
- fNavigators.push_back(navi);
- }
-
- if (navi != gGeoManager->GetCurrentNavigator()) {
- LOG(fatal) << "ca::tools::MaterialHelper: unexpected behavior of TGeoManager::AddNavigator() !!";
- }
- }
-
- if (!navi) {
- LOG(fatal) << "ca::tools::MaterialHelper: can not find / create TGeoNavigator for thread " << iThread;
- }
-
- return navi;
-}
-
-
-void MaterialHelper::CleanUpThreads()
-{
- // clean up thread data in TGeoManager
-
- gGeoManager->ClearThreadsMap();
-
- // create a default navigator for multithreaded mode
- // otherwise gGeoManager->GetCurrentNavigator() returns nullptr that might produces crashes in other code
-
- GetCurrentNavigator(-1);
-}
-
-void MaterialHelper::InitThreads()
-{
- // (re)initialise the number of threads in TGeoManager
-
- // reserve enough threads in TGeoManager
- if (fNthreads > gGeoManager->GetMaxThreads()) {
- gGeoManager->SetMaxThreads(fNthreads);
- // in case the number of threads is truncated by TGeoManager (must not happen)
- fNthreads = gGeoManager->GetMaxThreads();
- if (fNthreads < 1) {
- fNthreads = 1;
- }
- }
- CleanUpThreads();
-}
-
-
-MaterialMap MaterialHelper::GenerateMaterialMap(double zRef, double zMin, double zMax, double xyMax, int nBinsDim)
-{
- if (fDoRadialProjection) {
- if (fTargetZ >= zRef - 1.) {
- LOG(fatal) << "ca::tools::MaterialHelper: the material reference Z must be at least 1 cm downstream the target: "
- "target Z = "
- << fTargetZ << ", material reference z = " << zRef;
- }
- if (zMin < fTargetZ) {
- zMin = fTargetZ;
- }
- }
-
- if (!(zMin <= zRef && zRef <= zMax)) {
- LOG(fatal) << "ca::tools::MaterialHelper: the material parameters are inconsistent. It must be: zMin " << zMin
- << " <= zRef " << zRef << " <= zMax " << zMax;
- }
-
- MaterialMap matBudget;
- matBudget.Initialize(nBinsDim, xyMax, zRef, zMin, zMax);
- double binSizeX = 2. * xyMax / nBinsDim;
- double binSizeY = 2. * xyMax / nBinsDim;
-
- // call it every time. for the case the number of threads was reset in some other code
- InitThreads();
-
- long nThreadCrosses[fNthreads];
- long nThreadRays[fNthreads];
- long nThreadRaysExpected = nBinsDim * nBinsDim * fNraysBinPerDim * fNraysBinPerDim / fNthreads;
- for (int i = 0; i < fNthreads; i++) {
- nThreadCrosses[i] = 0;
- nThreadRays[i] = 0;
- }
-
- LOG(info) << "ca::tools::MaterialHelper: Generate material map using " << fNthreads << " threads..";
-
- auto naviThread = [&](int iThread) {
- // create a navigator for this thread
- TGeoNavigator* navi = GetCurrentNavigator(iThread);
-
- bool doPrint = false; // flag for debug printing
-
- for (int iBinX = iThread; iBinX < nBinsDim; iBinX += fNthreads) {
- for (int iBinY = 0; iBinY < nBinsDim; ++iBinY) {
- double radThick = 0.;
- int nRays = 0;
- double d = 1. / (fNraysBinPerDim);
- for (int iStepX = 0; iStepX < fNraysBinPerDim; iStepX++) {
- for (int iStepY = 0; iStepY < fNraysBinPerDim; iStepY++) {
-
- // ray position at zRef
- double rayX = -xyMax + (iBinX + d * (iStepX + 0.5)) * binSizeX;
- double rayY = -xyMax + (iBinY + d * (iStepY + 0.5)) * binSizeY;
-
- // ray slopes
- double tx = 0.;
- double ty = 0.;
- double t = 1.;
- if (fDoRadialProjection) {
- tx = rayX / (zRef - fTargetZ);
- ty = rayY / (zRef - fTargetZ);
- t = sqrt(1. + tx * tx + ty * ty);
- }
-
- // ray position at zMin
- double z = zMin;
- double x = rayX + tx * (z - zRef);
- double y = rayY + ty * (z - zRef);
- if (doPrint) {
- LOG(info) << "ray at " << x << " " << y << " " << z << " zMin " << zMin << " zMax " << zMax;
- }
-
- TGeoNode* node = navi->InitTrack(x, y, z, tx / t, ty / t, 1. / t);
-
- bool doContinue = 1;
- for (int iStep = 0; doContinue; iStep++) {
-
- if (!node) {
- // may happen when tracing outside of the CBM volume -> produce a warning
- LOG(warning) << "ca::tools::MaterialHelper: TGeoNavigator can not find the geo node";
- break;
- }
-
- TGeoMedium* medium = node->GetMedium();
- if (!medium) {
- LOG(fatal) << "ca::tools::MaterialHelper: TGeoNavigator can not find the geo medium";
- }
-
- TGeoMaterial* material = medium->GetMaterial();
- if (!material) {
- LOG(fatal) << "ca::tools::MaterialHelper: TGeoNavigator can not find the geo material";
- }
-
- double radLen = material->GetRadLen();
-
- if (radLen < 0.56) { // 0.5612 is rad. length of Lead at normal density
- LOG(fatal) << "ca::tools::MaterialHelper: failed assertion!";
- LOG(fatal) << " TGeoNavigator returns a suspicious material with an unrealistic "
- << "radiation length of " << radLen << " cm.";
- LOG(fatal) << " The allowed minimum is 0.56 cm (Lead has 0.5612 cm). Check your material! "
- "Modify this assertion if necessary.";
- LOG(fatal) << " TGeoNode \"" << node->GetName() << "\", TGeoMedium \"" << medium->GetName()
- << "\", TGeoMaterial \"" << material->GetName() << "\"";
- }
-
- // move navi to the next material border
- node = navi->FindNextBoundaryAndStep(); //5., kTRUE);
- nThreadCrosses[iThread]++;
-
- if (doPrint) {
- LOG(info) << " RL " << radLen << " next pt " << navi->GetCurrentPoint()[0] << " "
- << navi->GetCurrentPoint()[1] << " " << navi->GetCurrentPoint()[2];
- }
-
- double zNew = navi->GetCurrentPoint()[2];
- if (zNew >= zMax) {
- zNew = zMax;
- doContinue = 0;
- }
-
- //if (zNew < z) {
- //LOG(info) << " z curr " << z << " z new " << zNew << " dz " << zNew - z ; }
- double dz = zNew - z;
- if (dz < 0.) {
- if (iStep == 0 && dz > -1.e-8) {
- // when the ray starts exactly at the volume border, first dz might become negative,
- // probably due to some roundoff errors. So we let it be a bit negative for the first intersection.
- dz = 0.;
- }
- else {
- LOG(fatal)
- << "ca::tools::MaterialHelper: TGeoNavigator propagates the ray upstream. Something is wrong."
- << " z old " << z << " z new " << zNew << ", dz " << dz;
- }
- }
- radThick += dz / radLen;
- z = zNew;
- }
-
- nRays++;
- nThreadRays[iThread]++;
- if (nThreadRays[iThread] % 1000000 == 0) {
- LOG(info) << "ca::tools::MaterialHelper: report from thread " << iThread << ": material map is "
- << 100. * nThreadRays[iThread] / nThreadRaysExpected << " \% done";
- }
- }
- }
- radThick /= nRays;
- if (doPrint) {
- LOG(info) << " radThick " << radThick;
- }
- //doPrint = (radThick > 0.01);
- matBudget.SetRadThickBin(iBinX, iBinY, radThick);
- } // iBinX
- } // iBinY
- };
-
- std::vector<std::thread> threads(fNthreads);
-
- // run n threads
- for (int i = 0; i < fNthreads; i++) {
- threads[i] = std::thread(naviThread, i);
- }
-
- // wait for the threads to finish
- for (auto& th : threads) {
- th.join();
- }
-
- CleanUpThreads();
-
- long nRays = 0;
- long nCrosses = 0;
- for (int i = 0; i < fNthreads; i++) {
- nRays += nThreadRays[i];
- nCrosses += nThreadCrosses[i];
- }
-
- LOG(info) << "ca::tools::MaterialHelper: shooted " << nRays << " rays. Each ray crossed " << nCrosses * 1. / nRays
- << " material boundaries in average";
-
- return matBudget;
-}
diff --git a/reco/L1/catools/CaToolsMaterialHelper.h b/reco/L1/catools/CaToolsMaterialHelper.h
deleted file mode 100644
index 8a42ba95b76acb21b233083f923b3ca5c42e4263..0000000000000000000000000000000000000000
--- a/reco/L1/catools/CaToolsMaterialHelper.h
+++ /dev/null
@@ -1,94 +0,0 @@
-/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
- SPDX-License-Identifier: GPL-3.0-only
- Authors: Sergey Gorbunov [committer] */
-
-
-/// \file CaToolsMaterialHelper.h
-/// \brief A helper class to create material maps (header)
-/// \since 29.06.2023
-/// \author S.Gorbunov
-
-
-#ifndef CaToolsMaterialHelper_H
-#define CaToolsMaterialHelper_H
-
-#include "CaMaterialMap.h"
-#include "Rtypes.h"
-#include "TObject.h"
-
-class TGeoNavigator;
-
-namespace cbm::ca::tools
-{
- /// Class ca::tools::Debugger helps to debug CA tracking
- ///
- ///
- /// class to create ca::MaterialMap material maps form the ROOT geometry
- ///
- class MaterialHelper : public TObject {
- public:
- MaterialHelper();
- ~MaterialHelper();
-
- ///
- /// project rays radially from the targetZ througth the XY-bins at a reference z.
- ///
- void SetDoRadialProjection(double targetZ)
- {
- fDoRadialProjection = true;
- fTargetZ = targetZ;
- }
-
- ///
- /// project rays horisontally along the Z axis (default)
- ///
- void SetDoHorisontalProjection() { fDoRadialProjection = false; }
-
- ///
- /// shoot nRaysDim * nRaysDim rays for each bin in the map
- ///
- void SetNraysPerDim(int nRaysDim) { fNraysBinPerDim = (nRaysDim > 0) ? nRaysDim : 1; }
-
-
- /// generate a material map, collecting the material between [zMin, zMax]
- /// with radial rays from (0,0,targetZ) througth the XY-bins at z == zRef.
- ///
- /// It creates a map with [nBinsDim x nBinsDim] bins, of a size of [+-xyMax, +-xyMax]
- /// shooting fNraysBinPerDim x fNraysBinPerDim through each bin
- ///
- /// The calculated radiation thickness is a projection of the rad.thickness along the ray onto the Z axis.
- /// RadThick = sum of (dZ / radLength) over ray trajectory pieces
- ///
- /// When doRadialProjection==false the rays are shoot horisontally along the Z axis
- ///
- cbm::algo::ca::MaterialMap GenerateMaterialMap(double zRef, double zMin, double zMax, double xyMax, int nBinsDim);
-
- void SetSafeMaterialInitialization(bool val = true) { fDoSafeInitialization = val; }
-
- private:
- /// initialise the necessary amount of threads in TGeoManager
- void InitThreads();
-
- /// Clean up TGeoManager: threadIds, create a defailt navigator
- void CleanUpThreads();
-
- /// Get navigator for current thread, create it if it doesnt exist.
- /// Produce an error when anything goes wrong
- TGeoNavigator* GetCurrentNavigator(int iThread);
-
- private:
- int fNthreadsOld{0}; // number of threads in TGeoManager before the helper was created
- int fNthreads{0}; // number of threads
- bool fDoRadialProjection{false}; // project rays horisontally along the Z axis (special mode)
- int fNraysBinPerDim{3}; // shoot fNraysBinPerDim * fNraysBinPerDim rays in each map bin
- double fTargetZ{0.}; // z of the target for the radial projection
- std::vector<TGeoNavigator*> fNavigators{}; // list of created navigators
- bool fDoSafeInitialization{false}; // perform slow but safe initialization
- // to get around the crashes in TGeoVoxelFinder
- ClassDef(MaterialHelper, 0);
- };
-
-} // namespace cbm::ca::tools
-
-
-#endif