From b61a4c8082ba6d49adb0514fc0600d22fd6aab0c Mon Sep 17 00:00:00 2001 From: "s.zharko@gsi.de" <s.zharko@gsi.de> Date: Tue, 20 Feb 2024 22:33:59 +0100 Subject: [PATCH] Introduction of the CbmGenerateMaterialMaps task --- macro/L1/configs/CMakeLists.txt | 1 + macro/L1/configs/material_budget_example.yaml | 11 + macro/run/run_reco.C | 13 +- reco/L1/CMakeLists.txt | 3 +- reco/L1/CbmL1.cxx | 2 +- reco/L1/CbmL1.h | 5 +- reco/L1/CbmTrackingDetectorInterfaceInit.cxx | 15 +- reco/L1/CbmTrackingDetectorInterfaceInit.h | 3 + reco/L1/L1LinkDef.h | 1 + .../CbmGenerateMaterialMaps.cxx | 194 ++++++++++++++++++ .../CbmGenerateMaterialMaps.h | 132 ++++++++++++ 11 files changed, 373 insertions(+), 7 deletions(-) create mode 100644 macro/L1/configs/material_budget_example.yaml create mode 100644 reco/L1/OffLineInterface/CbmGenerateMaterialMaps.cxx create mode 100644 reco/L1/OffLineInterface/CbmGenerateMaterialMaps.h diff --git a/macro/L1/configs/CMakeLists.txt b/macro/L1/configs/CMakeLists.txt index 3d65acd46c..79d7eea92c 100644 --- a/macro/L1/configs/CMakeLists.txt +++ b/macro/L1/configs/CMakeLists.txt @@ -3,5 +3,6 @@ Install(FILES config_ideal_hits_mcbm.yaml ca_params_sts.yaml ca_params_global.yaml ca_params_user_example.yaml + material_budget_example.yaml DESTINATION share/cbmroot/macro/L1/configs ) diff --git a/macro/L1/configs/material_budget_example.yaml b/macro/L1/configs/material_budget_example.yaml new file mode 100644 index 0000000000..8030bfa89d --- /dev/null +++ b/macro/L1/configs/material_budget_example.yaml @@ -0,0 +1,11 @@ +--- +user_slices: + - { ref_z: 10, min_z: 0, max_z: 20, max_xy: 100, name: "test_0" } + - { ref_z: 30, min_z: 20, max_z: 40, max_xy: 100, name: "test_1" } +pitch: 0.05 +max_nof_bins: 100 +nof_rays: 3 +parallel_rays: true +tracking_stations: true +safe_material_init: false +... diff --git a/macro/run/run_reco.C b/macro/run/run_reco.C index 455dfa21e9..4fc60cc382 100644 --- a/macro/run/run_reco.C +++ b/macro/run/run_reco.C @@ -17,6 +17,7 @@ #include "CbmDefs.h" #include "CbmFindPrimaryVertex.h" #include "CbmFsdHitProducer.h" +#include "CbmGenerateMaterialMaps.h" #include "CbmKF.h" #include "CbmL1.h" #include "CbmL1StsTrackFinder.h" @@ -412,15 +413,25 @@ void run_reco(TString input = "", Int_t nTimeSlices = -1, Int_t firstTimeSlice = // ----- Track finding in STS (+ MVD) -------------------------------- if (useMvd || useSts) { - auto pDetIF = new CbmTrackingDetectorInterfaceInit(); + auto* pDetIF = new CbmTrackingDetectorInterfaceInit(); run->AddTask(pDetIF); // Geometry interface initializer for tracker std::cout << "-I- " << myName << ": Added task " << pDetIF->GetName() << std::endl; + // Kalman filter auto kalman = new CbmKF(); run->AddTask(kalman); std::cout << "-I- " << myName << ": Added task " << kalman->GetName() << std::endl; + // Material map generator + { + TString cfgFile = srcDir + "/macro/L1/configs/material_budget_example.yaml"; + auto* pMaterialMapGenerator = new CbmGenerateMaterialMaps(); + pMaterialMapGenerator->SetUserConfig(cfgFile); + pMaterialMapGenerator->SetOutputName("material.root"); + //run->AddTask(pMaterialMapGenerator); + } + if (debugWithMC) { auto* pIdealHitProducer = new cbm::ca::IdealHitProducer("IdealHitProducer", 3); pIdealHitProducer->SetConfigName(srcDir + "/macro/L1/configs/config_ideal_hits.yaml"); diff --git a/reco/L1/CMakeLists.txt b/reco/L1/CMakeLists.txt index 14f11c0ee2..8ea34bed88 100644 --- a/reco/L1/CMakeLists.txt +++ b/reco/L1/CMakeLists.txt @@ -31,6 +31,7 @@ set(SRCS OffLineInterface/CbmL1StsTrackFinder.cxx OffLineInterface/CbmL1GlobalTrackFinder.cxx OffLineInterface/CbmL1GlobalFindTracksEvents.cxx + OffLineInterface/CbmGenerateMaterialMaps.cxx CbmL1Util.cxx CbmL1Performance.cxx @@ -51,7 +52,7 @@ set(SRCS catools/CaToolsWindowFinder.cxx catools/CaToolsWFExpression.cxx catools/CaToolsMaterialHelper.cxx - + qa/CbmTrackerInputQaTrd.cxx qa/CbmTrackerInputQaTof.cxx qa/CbmCaInputQaBase.cxx diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx index 4696aea698..5c66cfaa0e 100644 --- a/reco/L1/CbmL1.cxx +++ b/reco/L1/CbmL1.cxx @@ -816,7 +816,7 @@ void CbmL1::GenerateMaterialMaps() } double zNew = 0.5 * (z1 + z2); - double maxXY = 1.3 * std::max(pStation->GetXmax(), pStation->GetYmax()); + double maxXY = 1.3 * std::max(std::fabs(pStation->GetXmax()), std::fabs(pStation->GetYmax())); //double maxXY = 80; // calculate n bins from the minimal pitch int nBins = static_cast<int>(std::ceil(2. * maxXY / fMatBudgetPitch)); diff --git a/reco/L1/CbmL1.h b/reco/L1/CbmL1.h index fb798ece09..67a7be9285 100644 --- a/reco/L1/CbmL1.h +++ b/reco/L1/CbmL1.h @@ -298,6 +298,9 @@ class CbmL1 : public FairTask { << " ! The material budget files are not used anymore !"; } + /// A helper for GetTargetInfo() + static void FindTargetNode(TString& targetPath, TGeoNode*& targetNode); + private: struct TH1FParameters { TString name, title; @@ -395,8 +398,6 @@ class CbmL1 : public FairTask { /// Get the target information void GetTargetInfo(); - /// A helper for GetTargetInfo() - void FindTargetNode(TString& targetPath, TGeoNode*& targetNode); void WriteHistosCurFile(TObject* obj); diff --git a/reco/L1/CbmTrackingDetectorInterfaceInit.cxx b/reco/L1/CbmTrackingDetectorInterfaceInit.cxx index 13a88fa821..f76300db64 100644 --- a/reco/L1/CbmTrackingDetectorInterfaceInit.cxx +++ b/reco/L1/CbmTrackingDetectorInterfaceInit.cxx @@ -22,7 +22,7 @@ ClassImp(CbmTrackingDetectorInterfaceInit) CbmTrackingDetectorInterfaceInit* CbmTrackingDetectorInterfaceInit::fpInstance = nullptr; -//------------------------------------------------------------------------------------------------------------------------------------- +// --------------------------------------------------------------------------------------------------------------------- // CbmTrackingDetectorInterfaceInit::CbmTrackingDetectorInterfaceInit() : FairTask("CbmTrackingDetectorInterfaceInit") { @@ -72,7 +72,7 @@ CbmTrackingDetectorInterfaceInit::CbmTrackingDetectorInterfaceInit() : FairTask( } } -//------------------------------------------------------------------------------------------------------------------------------------- +// --------------------------------------------------------------------------------------------------------------------- // CbmTrackingDetectorInterfaceInit::~CbmTrackingDetectorInterfaceInit() { @@ -80,3 +80,14 @@ CbmTrackingDetectorInterfaceInit::~CbmTrackingDetectorInterfaceInit() fpInstance = nullptr; } } + +// --------------------------------------------------------------------------------------------------------------------- +// +std::vector<const CbmTrackingDetectorInterfaceBase*> CbmTrackingDetectorInterfaceInit::GetActiveInterfaces() const +{ + std::vector<const CbmTrackingDetectorInterfaceBase*> vRes; + for (const TObject* pTask : *(this->GetListOfTasks())) { + vRes.push_back(dynamic_cast<const CbmTrackingDetectorInterfaceBase*>(pTask)); + } + return vRes; +} diff --git a/reco/L1/CbmTrackingDetectorInterfaceInit.h b/reco/L1/CbmTrackingDetectorInterfaceInit.h index 7cb01e3248..f465cbc394 100644 --- a/reco/L1/CbmTrackingDetectorInterfaceInit.h +++ b/reco/L1/CbmTrackingDetectorInterfaceInit.h @@ -20,6 +20,7 @@ class CbmStsTrackingInterface; class CbmMuchTrackingInterface; class CbmTrdTrackingInterface; class CbmTofTrackingInterface; +class CbmTrackingDetectorInterfaceBase; /// Class CbmTrackingDetectorInterfaceInit implements a task for the tracking detector interfaces initialization. /// The tracking detector interfaces are added as subtasks. The CbmTrackingDetectorInterfaceInit class instance is to @@ -50,6 +51,8 @@ class CbmTrackingDetectorInterfaceInit : public FairTask { CbmTrdTrackingInterface* GetTrdTrackingInterface() { return fpTrdTrackingInterface; } CbmTofTrackingInterface* GetTofTrackingInterface() { return fpTofTrackingInterface; } + /// \brief Returns a vector of pointers to the tracking detector interfaces, included into the setup + std::vector<const CbmTrackingDetectorInterfaceBase*> GetActiveInterfaces() const; private: static CbmTrackingDetectorInterfaceInit* fpInstance; ///< Instance of the class diff --git a/reco/L1/L1LinkDef.h b/reco/L1/L1LinkDef.h index 34a671499a..84bd9fd6de 100644 --- a/reco/L1/L1LinkDef.h +++ b/reco/L1/L1LinkDef.h @@ -19,6 +19,7 @@ //#pragma link C++ class CbmL1MuchFinderQa+; #pragma link C++ class CbmL1GlobalTrackFinder + ; #pragma link C++ class CbmL1GlobalFindTracksEvents + ; +#pragma link C++ class CbmGenerateMaterialMaps.h + ; //#pragma link C++ class CbmL1SttHit+; //#pragma link C++ class CbmL1SttTrackFinder+; //#pragma link C++ class CbmL1SttTrack+; diff --git a/reco/L1/OffLineInterface/CbmGenerateMaterialMaps.cxx b/reco/L1/OffLineInterface/CbmGenerateMaterialMaps.cxx new file mode 100644 index 0000000000..07e57c83d1 --- /dev/null +++ b/reco/L1/OffLineInterface/CbmGenerateMaterialMaps.cxx @@ -0,0 +1,194 @@ +/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt + SPDX-License-Identifier: GPL-3.0-only + Authors: Sergei Zharko [committer] */ + +/// \file CbmGenerateMaterialMaps.h +/// \brief A FAIR-task to generate material maps (header) +/// \since 18.02.2024 +/// \author Sergei Zharko <s.zharko@gsi.de> + +#include "CbmGenerateMaterialMaps.h" + +#include "CbmL1.h" // for FindTargetNode +#include "CbmTrackingDetectorInterfaceBase.h" +#include "CbmTrackingDetectorInterfaceInit.h" +#include "TFile.h" +#include "TH2F.h" + +// --------------------------------------------------------------------------------------------------------------------- +// +CbmGenerateMaterialMaps::CbmGenerateMaterialMaps() : CbmGenerateMaterialMaps("CbmGenerateMaterialMaps", 0) {} + +// --------------------------------------------------------------------------------------------------------------------- +// +CbmGenerateMaterialMaps::CbmGenerateMaterialMaps(const char* name, int verbose) : FairTask(name, verbose) {} + +// --------------------------------------------------------------------------------------------------------------------- +// +void CbmGenerateMaterialMaps::Finish() {} + +// --------------------------------------------------------------------------------------------------------------------- +// +InitStatus CbmGenerateMaterialMaps::Init() +{ + // ----- Init configuration + if (!fsUserConfig.empty()) { + fConfig = cbm::algo::config::ReadFromFile<Config>(fsUserConfig); + } + else { + LOG(info) << fName << ": configuration file was not provided. Using default settings"; + } + + InitTarget(); + + fpMaterialHelper = std::make_unique<cbm::ca::tools::MaterialHelper>(); + if (!fConfig.fbParallelRays) { + fpMaterialHelper->SetDoRadialProjection(fTargetZ); + } + fpMaterialHelper->SetSafeMaterialInitialization(true); + fpMaterialHelper->SetNraysPerDim(fConfig.fNofRays); + + if (fConfig.fbTrackingStations) { + std::set<MaterialSlice> mSlice; + // Loop over stations + auto* pIfs = CbmTrackingDetectorInterfaceInit::Instance(); + if (pIfs) { + for (const auto* pIf : pIfs->GetActiveInterfaces()) { + auto detName = pIf->GetDetectorName(); + LOG(info) << fName << ": Generating material budget map for " << detName; + int nSt = pIf->GetNtrackingStations(); + for (int iSt = 0; iSt < nSt; ++iSt) { + MaterialSlice slice; + slice.fName = detName + Form("_station_%d", iSt); + slice.fRefZ = pIf->GetZref(iSt); + slice.fMinZ = pIf->GetZmin(iSt); + slice.fMaxZ = pIf->GetZmax(iSt); + slice.fMaxXY = std::max(std::fabs(pIf->GetXmax(iSt)), std::fabs(pIf->GetYmax(iSt))); + mSlice.insert(slice); + } + } + double zLast = fTargetZ + 1.; + // Run material helper + for (auto it = mSlice.begin(); it != mSlice.end(); ++it) { + LOG(info) << "Creating material map for " << it->fName; + double z1 = it->fMaxZ; + double z2 = z1; + auto itNext = std::next(it, 1); + if (itNext != mSlice.end()) { + z2 = itNext->fMinZ; + } + double zRef = it->fRefZ; + double zNew = 0.5 * (z1 + z2); + double xyMax = kXYoffset * it->fMaxXY; + int nBins = static_cast<int>(std::ceil(2. * xyMax / fConfig.fPitch)); + if (nBins < 1) { + LOG(fatal) << fName << ": selected pitch " << fConfig.fPitch << " gives wrong number of bins = " << nBins; + } + if (nBins > fConfig.fMaxNofBins) { + nBins = fConfig.fMaxNofBins; + } + fmMaterial[it->fName] = fpMaterialHelper->GenerateMaterialMap(zRef, zLast, zNew, xyMax, nBins); + zLast = zNew; + } + } + else { + LOG(error) << fName << ": tracking detector interfaces are not defined, so the materials maps cannot be " + << "generated for the tracking stations"; + } + } + + // Loop over user-defined slices + for (const auto& slice : fConfig.fvUserSlices) { + LOG(info) << "Creating material map for " << slice.fName; + if (fmMaterial.find(slice.fName) != fmMaterial.end()) { + LOG(warn) << fName << ": Material for slice " << slice.fName << " was already prepared. " + << "Please, check your configuration file"; + continue; + } + double xyMax = kXYoffset * slice.fMaxXY; + double zMin = slice.fMinZ; + double zMax = slice.fMaxZ; + double zRef = slice.fRefZ; + + int nBins = static_cast<int>(std::ceil(2. * xyMax / fConfig.fPitch)); + + if (nBins < 1) { + LOG(fatal) << fName << ": selected pitch " << fConfig.fPitch << " gives wrong number of bins = " << nBins; + } + if (nBins > fConfig.fMaxNofBins) { + nBins = fConfig.fMaxNofBins; + } + fmMaterial[slice.fName] = fpMaterialHelper->GenerateMaterialMap(zRef, zMin, zMax, xyMax, nBins); + } + + WriteMaterialMaps(); + return kSUCCESS; +} + +// --------------------------------------------------------------------------------------------------------------------- +// +InitStatus CbmGenerateMaterialMaps::ReInit() { return kSUCCESS; } + +// --------------------------------------------------------------------------------------------------------------------- +// +void CbmGenerateMaterialMaps::InitTarget() +{ + // Loop over all nodes till a node with name "target" is found + // Extract the required infrmation from the node and store it in the + // proper structure + // The complete logic depends on the naming convention of the target. + // If the node doesn't contain the string target the procedure will fail + + fTargetX = 1.e10; + fTargetY = 1.e10; + fTargetZ = 1.e10; + + TString targetPath; + TGeoNode* targetNode{nullptr}; + CbmL1::FindTargetNode(targetPath, targetNode); + + if (!targetNode) { + LOG(fatal) << fName << ": target node is not found"; + } + + 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); + fTargetX = global[0]; + fTargetY = global[1]; + fTargetZ = global[2]; + + LOG(info) << fName << ": target found: \"" << targetPath << "\" at ( " << fTargetX << " " << fTargetY << " " + << fTargetZ << " ) [cm] "; +} + +// --------------------------------------------------------------------------------------------------------------------- +// +void CbmGenerateMaterialMaps::WriteMaterialMaps() +{ + auto f = TFile(fsOutputFile.c_str(), "RECREATE"); + f.cd(); + for (const auto& [name, material] : fmMaterial) { + double zMin = material.GetZmin(); + double zMax = material.GetZmax(); + std::string title = Form("Material thickness in X_{0}, z region [%.2f,%.2f] cm (%s)", zMin, zMax, name.c_str()); + if (fConfig.fbParallelRays) { + title += "; horizontal projection"; + } + else { + title += Form("; radial projection from z=%f cm", fTargetZ); + } + double nBins = material.GetNbins(); + double xyMax = material.GetXYmax(); + auto* h = new TH2F(name.c_str(), title.c_str(), nBins, -xyMax, xyMax, nBins, -xyMax, xyMax); + h->GetXaxis()->SetTitle("x [cm]"); + h->GetYaxis()->SetTitle("y [cm]"); + for (int iBinX = 0; iBinX < nBins; ++iBinX) { + for (int iBinY = 0; iBinY < nBins; ++iBinY) { + h->SetBinContent(iBinX + 1, iBinY + 1, material.GetRadThickBin(iBinX, iBinY)); + } + } + } + f.Write(); +} diff --git a/reco/L1/OffLineInterface/CbmGenerateMaterialMaps.h b/reco/L1/OffLineInterface/CbmGenerateMaterialMaps.h new file mode 100644 index 0000000000..bde2425c7f --- /dev/null +++ b/reco/L1/OffLineInterface/CbmGenerateMaterialMaps.h @@ -0,0 +1,132 @@ +/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt + SPDX-License-Identifier: GPL-3.0-only + Authors: Sergei Zharko [committer] */ + +/// \file CbmGenerateMaterialMaps.h +/// \brief A FAIR-task to generate material maps (header) +/// \since 18.02.2024 +/// \author Sergei Zharko <s.zharko@gsi.de> + +#pragma once + +#include "CaMaterialMap.h" +#include "CaToolsMaterialHelper.h" +#include "FairTask.h" +#include "config/Yaml.h" + +#include <map> +#include <string> +#include <string_view> +#include <tuple> +#include <vector> + +/// \class CbmGenerateMaterialMaps +/// \brief Steer class for executing the material budget maps generator independently from tracking +class CbmGenerateMaterialMaps : public FairTask { + using MaterialMap = cbm::algo::ca::MaterialMap; + using MaterialHelper = cbm::ca::tools::MaterialHelper; + + /// \struct MaterialSlice + /// \brief Input parameters for the material map generation + struct MaterialSlice { + std::string fName; + double fRefZ; ///< Reference z-coordinate [cm] + double fMinZ; ///< Lower bound along z-axis [cm] + double fMaxZ; ///< Upper bound along z-axis [cm] + double fMaxXY; ///< Size in the transverse plane [cm] + + bool operator<(const MaterialSlice& r) const { return (fRefZ < r.fRefZ); } + + CBM_PROPERTIES(cbm::algo::config::Property(&CbmGenerateMaterialMaps::MaterialSlice::fName, "name", ""), + cbm::algo::config::Property(&CbmGenerateMaterialMaps::MaterialSlice::fRefZ, "ref_z", ""), + cbm::algo::config::Property(&CbmGenerateMaterialMaps::MaterialSlice::fMinZ, "min_z", ""), + cbm::algo::config::Property(&CbmGenerateMaterialMaps::MaterialSlice::fMaxZ, "max_z", ""), + cbm::algo::config::Property(&CbmGenerateMaterialMaps::MaterialSlice::fMaxXY, "max_xy", "")); + }; + + struct Config { + std::vector<MaterialSlice> fvUserSlices; ///< Material slices defined by user + + double fPitch = 0.1; ///< Minimal bin size (cm) + int fMaxNofBins = 100; ///< Number of bins in material budged map (x and y axes) + int fNofRays = 3; ///< Number of rays per dimension in each bin + bool fbParallelRays = false; ///< Rays mode (false - radial, true - parallel to z-axis) + bool fbTrackingStations = false; ///< Generates material maps for the actual geometry stations + bool fbSafeMaterialInit = true; ///< Safe material initialization (takes extra computational time) + + CBM_PROPERTIES(cbm::algo::config::Property(&CbmGenerateMaterialMaps::Config::fvUserSlices, "user_slices", ""), + cbm::algo::config::Property(&CbmGenerateMaterialMaps::Config::fPitch, "pitch", ""), + cbm::algo::config::Property(&CbmGenerateMaterialMaps::Config::fMaxNofBins, "max_nof_bins", ""), + cbm::algo::config::Property(&CbmGenerateMaterialMaps::Config::fNofRays, "nof_rays", ""), + cbm::algo::config::Property(&CbmGenerateMaterialMaps::Config::fbParallelRays, "parallel_rays", ""), + cbm::algo::config::Property(&CbmGenerateMaterialMaps::Config::fbTrackingStations, + "tracking_stations", ""), + cbm::algo::config::Property(&CbmGenerateMaterialMaps::Config::fbSafeMaterialInit, + "safe_material_init", "")); + }; + + public: + /// \brief Default constructor + CbmGenerateMaterialMaps(); + + /// \brief Constructor from parameters + /// \param name Name of the task + /// \param verbose Verbosity level + CbmGenerateMaterialMaps(const char* name, int verbose); + + /// \brief Destructor + virtual ~CbmGenerateMaterialMaps() = default; + + /// \brief Copy constructor + CbmGenerateMaterialMaps(const CbmGenerateMaterialMaps&) = delete; + + /// \brief Move constructor + CbmGenerateMaterialMaps(CbmGenerateMaterialMaps&&) = delete; + + /// \brief Copy assignment operator + CbmGenerateMaterialMaps& operator=(const CbmGenerateMaterialMaps&) = delete; + + /// \brief Move assignment operator + CbmGenerateMaterialMaps& operator=(CbmGenerateMaterialMaps&&) = delete; + + /// \brief Initializer + InitStatus Init(); + + /// \brief Re-initializer + InitStatus ReInit(); + + /// \brief Finish function + void Finish(); + + /// \brief Sets output file name + void SetOutputName(const char* fileName) { fsOutputFile = fileName; } + + /// \brief Sets user configuration file + void SetUserConfig(const char* userConfig) { fsUserConfig = userConfig; } + + /// \brief Gets material map + const std::map<std::string, MaterialMap>& GetMaterial() const { return fmMaterial; } + + private: + /// \brief Gets target + void InitTarget(); + + /// \brief Writes material budget maps to file + /// + /// The material maps are represented with ROOT histograms. The value of the material thickness + /// is represented in units of radiational length + void WriteMaterialMaps(); + + static constexpr double kXYoffset = 1.3; /// Offset from station boarders [scale] + static constexpr double kTargetOffset = 1.; /// Offset from target to start generating mat. maps [cm] + + Config fConfig; ///< Configuration for the task + std::string fsUserConfig = ""; ///< User configuration file + std::string fsOutputFile = "matBudget.root"; ///< Output file name + std::map<std::string, MaterialMap> fmMaterial; ///< Material budget maps + std::unique_ptr<MaterialHelper> fpMaterialHelper = nullptr; ///< Material helper instance + + double fTargetX; + double fTargetY; + double fTargetZ; +}; -- GitLab