From 94d20193a181e658c0938d8ff6a0226790ee4a6d Mon Sep 17 00:00:00 2001 From: "se.gorbunov" <se.gorbunov@gsi.de> Date: Fri, 30 Jun 2023 00:31:29 +0000 Subject: [PATCH] L1: get rid of material budget map files. Generate maps on-the-fly when initialising the tracker --- reco/L1/CMakeLists.txt | 4 + reco/L1/CbmL1.cxx | 204 ++++++++++----- reco/L1/CbmL1.h | 30 ++- reco/L1/L1Algo/L1BaseStationInfo.cxx | 2 - reco/L1/L1Algo/L1Material.cxx | 156 +++-------- reco/L1/L1Algo/L1Material.h | 52 +++- reco/L1/L1Algo/L1MaterialMonitor.cxx | 146 +++++++++++ reco/L1/L1Algo/L1MaterialMonitor.h | 105 ++++++++ reco/L1/L1LinkDef.h | 1 + reco/L1/catools/CaToolsMaterialHelper.cxx | 304 ++++++++++++++++++++++ reco/L1/catools/CaToolsMaterialHelper.h | 95 +++++++ 11 files changed, 892 insertions(+), 207 deletions(-) create mode 100644 reco/L1/L1Algo/L1MaterialMonitor.cxx create mode 100644 reco/L1/L1Algo/L1MaterialMonitor.h create mode 100644 reco/L1/catools/CaToolsMaterialHelper.cxx create mode 100644 reco/L1/catools/CaToolsMaterialHelper.h diff --git a/reco/L1/CMakeLists.txt b/reco/L1/CMakeLists.txt index af613208b7..48e538c94f 100644 --- a/reco/L1/CMakeLists.txt +++ b/reco/L1/CMakeLists.txt @@ -57,6 +57,7 @@ set(SRCS CbmL1MCTrack.cxx CbmL1Track.cxx L1Algo/L1Material.cxx + L1Algo/L1MaterialMonitor.cxx L1Algo/L1UMeasurementInfo.cxx L1Algo/L1XYMeasurementInfo.cxx L1Algo/L1Field.cxx @@ -79,6 +80,8 @@ set(SRCS catools/CaToolsDebugger.cxx catools/CaToolsWindowFinder.cxx catools/CaToolsWFExpression.cxx + catools/CaToolsMaterialHelper.cxx + ParticleFinder/CbmL1PFFitter.cxx ParticleFinder/CbmL1PFMCParticle.cxx @@ -120,6 +123,7 @@ set(HEADERS catools/CaToolsLinkKey.h catools/CaToolsHitRecord.h #utils/CbmCaIdealHitProducer.h + catools/CaToolsMaterialHelper.h qa/CbmCaInputQaBase.h ) diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx index 8b43f3d064..c9db466cb9 100644 --- a/reco/L1/CbmL1.cxx +++ b/reco/L1/CbmL1.cxx @@ -32,6 +32,8 @@ #include "CbmTrdTrackingInterface.h" #include <boost/filesystem.hpp> + +#include "CaToolsMaterialHelper.h" // TODO: include of CbmSetup.h creates problems on Mac // #include "CbmSetup.h" #include "CbmEvent.h" @@ -58,6 +60,7 @@ #include "TVectorD.h" #include <TFile.h> +#include <chrono> #include <fstream> #include <iomanip> #include <iostream> @@ -230,9 +233,9 @@ InitStatus CbmL1::Init() if (L1Algo::TrackingMode::kGlobal == fTrackingMode) { //at the moment trd2d tracking only - fUseMVD = 0; - fUseSTS = 0; - fUseMUCH = 0; + fUseMVD = 1; + fUseSTS = 1; + fUseMUCH = 1; fUseTRD = 1; fUseTOF = 0; @@ -464,9 +467,11 @@ InitStatus CbmL1::Init() // ** Stations geometry initialization ** // *************************************** + std::vector<L1BaseStationInfo> vStations; + vStations.reserve(100); + // *** MVD stations info *** if (fUseMVD) { - auto materialTableMvd = ReadMaterialBudget(L1DetectorID::kMvd); for (int iSt = 0; iSt < fNMvdStationsGeom; ++iSt) { auto stationInfo = L1BaseStationInfo(L1DetectorID::kMvd, iSt); stationInfo.SetStationType(1); // MVD @@ -478,7 +483,6 @@ InitStatus CbmL1::Init() stationInfo.SetRmin(mvdInterface->GetRmin(iSt)); stationInfo.SetRmax(mvdInterface->GetRmax(iSt)); stationInfo.SetZthickness(mvdInterface->GetThickness(iSt)); - stationInfo.SetMaterialMap(std::move(materialTableMvd[iSt])); // TODO: The CA TF result is dependent from type of geometry settings. Should be understood (S.Zharko) stationInfo.SetFrontBackStripsGeometry((fscal) mvdInterface->GetStripsStereoAngleFront(iSt), (fscal) mvdInterface->GetStripsStereoAngleBack(iSt)); @@ -486,14 +490,13 @@ InitStatus CbmL1::Init() if (fvmDisabledStationIDs[L1DetectorID::kMvd].find(iSt) != fvmDisabledStationIDs[L1DetectorID::kMvd].cend()) { stationInfo.SetTrackingStatus(false); } - fInitManager.AddStation(stationInfo); - LOG(info) << "- MVD station " << iSt << " at z = " << stationInfo.GetZdouble() << " cm"; + LOG(info) << "- MVD station " << iSt << " at z = " << stationInfo.GetZdouble() << " cm "; + vStations.push_back(stationInfo); } } // *** STS stations info *** if (fUseSTS) { - auto materialTableSts = ReadMaterialBudget(L1DetectorID::kSts); for (int iSt = 0; iSt < fNStsStationsGeom; ++iSt) { auto stationInfo = L1BaseStationInfo(L1DetectorID::kSts, iSt); stationInfo.SetStationType(0); // STS @@ -507,7 +510,7 @@ InitStatus CbmL1::Init() stationInfo.SetRmin(stsInterface->GetRmin(iSt)); stationInfo.SetRmax(stsInterface->GetRmax(iSt)); stationInfo.SetZthickness(stsInterface->GetThickness(iSt)); - stationInfo.SetMaterialMap(std::move(materialTableSts[iSt])); + // TODO: The CA TF result is dependent from type of geometry settings. Should be understood (S.Zharko) stationInfo.SetFrontBackStripsGeometry((fscal) stsInterface->GetStripsStereoAngleFront(iSt), (fscal) stsInterface->GetStripsStereoAngleBack(iSt)); @@ -515,14 +518,13 @@ InitStatus CbmL1::Init() if (fvmDisabledStationIDs[L1DetectorID::kSts].find(iSt) != fvmDisabledStationIDs[L1DetectorID::kSts].cend()) { stationInfo.SetTrackingStatus(false); } - fInitManager.AddStation(stationInfo); - LOG(info) << "- STS station " << iSt << " at z = " << stationInfo.GetZdouble() << " cm"; + LOG(info) << "- STS station " << iSt << " at z = " << stationInfo.GetZdouble() << " cm "; + vStations.push_back(stationInfo); } } // *** MuCh stations info *** if (fUseMUCH) { - auto materialTableMuch = ReadMaterialBudget(L1DetectorID::kMuch); for (int iSt = 0; iSt < fNMuchStationsGeom; ++iSt) { auto stationInfo = L1BaseStationInfo(L1DetectorID::kMuch, iSt); stationInfo.SetStationType(2); // MuCh @@ -536,7 +538,6 @@ InitStatus CbmL1::Init() stationInfo.SetRmin(muchInterface->GetRmin(iSt)); stationInfo.SetRmax(muchInterface->GetRmax(iSt)); stationInfo.SetZthickness(muchInterface->GetThickness(iSt)); - stationInfo.SetMaterialMap(std::move(materialTableMuch[iSt])); // TODO: The CA TF result is dependent from type of geometry settings. Should be understood (S.Zharko) stationInfo.SetFrontBackStripsGeometry((fscal) muchInterface->GetStripsStereoAngleFront(iSt), (fscal) muchInterface->GetStripsStereoAngleBack(iSt)); @@ -544,14 +545,13 @@ InitStatus CbmL1::Init() if (fvmDisabledStationIDs[L1DetectorID::kMuch].find(iSt) != fvmDisabledStationIDs[L1DetectorID::kMuch].cend()) { stationInfo.SetTrackingStatus(false); } - fInitManager.AddStation(stationInfo); + vStations.push_back(stationInfo); LOG(info) << "- MuCh station " << iSt << " at z = " << stationInfo.GetZdouble() << " cm"; } } // *** TRD stations info *** if (fUseTRD) { - auto materialTableTrd = ReadMaterialBudget(L1DetectorID::kTrd); for (int iSt = 0; iSt < fNTrdStationsGeom; ++iSt) { auto stationInfo = L1BaseStationInfo(L1DetectorID::kTrd, iSt); stationInfo.SetStationType((iSt == 1 || iSt == 3) ? 6 : 3); // MuCh @@ -565,7 +565,6 @@ InitStatus CbmL1::Init() stationInfo.SetRmin(trdInterface->GetRmin(iSt)); stationInfo.SetRmax(trdInterface->GetRmax(iSt)); stationInfo.SetZthickness(trdInterface->GetThickness(iSt)); - stationInfo.SetMaterialMap(std::move(materialTableTrd[iSt])); fscal trdFrontPhi = trdInterface->GetStripsStereoAngleFront(iSt); fscal trdBackPhi = trdInterface->GetStripsStereoAngleBack(iSt); if (L1Algo::TrackingMode::kGlobal == fTrackingMode) { stationInfo.SetTimeInfo(false); } @@ -577,14 +576,13 @@ InitStatus CbmL1::Init() if (fvmDisabledStationIDs[L1DetectorID::kTrd].find(iSt) != fvmDisabledStationIDs[L1DetectorID::kTrd].cend()) { stationInfo.SetTrackingStatus(false); } - fInitManager.AddStation(stationInfo); + vStations.push_back(stationInfo); LOG(info) << "- TRD station " << iSt << " at z = " << stationInfo.GetZdouble() << " cm"; } } // *** TOF stations info *** if (fUseTOF) { - auto materialTableTof = ReadMaterialBudget(L1DetectorID::kTof); for (int iSt = 0; iSt < fNTofStationsGeom; ++iSt) { auto stationInfo = L1BaseStationInfo(L1DetectorID::kTof, iSt); stationInfo.SetStationType(4); @@ -595,7 +593,6 @@ InitStatus CbmL1::Init() stationInfo.SetZ(tofInterface->GetZ(iSt)); auto thickness = tofInterface->GetThickness(iSt); stationInfo.SetZthickness(thickness); - stationInfo.SetMaterialMap(std::move(materialTableTof[iSt])); stationInfo.SetXmax(tofInterface->GetXmax(iSt)); stationInfo.SetYmax(tofInterface->GetYmax(iSt)); stationInfo.SetRmin(tofInterface->GetRmin(iSt)); @@ -607,11 +604,63 @@ InitStatus CbmL1::Init() if (fvmDisabledStationIDs[L1DetectorID::kTof].find(iSt) != fvmDisabledStationIDs[L1DetectorID::kTof].cend()) { stationInfo.SetTrackingStatus(false); } - fInitManager.AddStation(stationInfo); + vStations.push_back(stationInfo); LOG(info) << "- TOF station " << iSt << " at z = " << stationInfo.GetZdouble() << " cm"; } } + // **************************************** + // ** ** + // ** Initialize material maps ** + // ** ** + // **************************************** + + { + LOG(info) << "- Generate material maps.."; + auto timerStart = std::chrono::high_resolution_clock::now(); + + ca::tools::MaterialHelper matHelper; + if (!fMatBudgetParallelProjection) { matHelper.SetDoRadialProjection(target.z); } + matHelper.SetNraysPerDim(fMatBudgetNrays); + + double zLast = target.z + 1.; // some gap (+1cm) to skip the target material + + for (unsigned int ist = 0; ist < vStations.size(); ist++) { + auto& station = vStations[ist]; + double z1 = station.GetZdouble() + station.GetZthickness()[0] / 2.; + double z2 = z1; + if (ist < vStations.size() - 1) { + // split materials between the stations at the middle + auto& stationNext = vStations[ist + 1]; + z2 = stationNext.GetZdouble() - stationNext.GetZthickness()[0] / 2.; + } + double zNew = 0.5 * (z1 + z2); + + double maxXY = 1.3 * station.GetRmax()[0]; + //double maxXY = 80; + // calculate n bins from the minimal pitch + int nBins = static_cast<int>(std::ceil(2. * maxXY / fMatBudgetPitch)); + if (nBins < 1) { LOG(fatal) << " material nBins " << nBins << " is not positive, something is wrong"; } + if (nBins > fMatBudgetNbins) { nBins = fMatBudgetNbins; } + + L1Material matBudget = matHelper.GenerateMaterialMap(station.GetZdouble(), zLast, zNew, maxXY, nBins); + + station.SetMaterialMap(matBudget); + + LOG(info) << "- Generated material map for tracking station " << ist << " at z = " << station.GetZdouble() + << " cm." + << " Material is collected between z = " << zLast << " and z = " << zNew; + + fInitManager.AddStation(station); + + zLast = zNew; + } + + auto timerEnd = std::chrono::high_resolution_clock::now(); + double materialGenerationTime = (double) (std::chrono::duration<double>(timerEnd - timerStart).count()); + LOG(info) << "- Generating material maps tooks " << materialGenerationTime << " seconds"; + } + // **************************************** // ** ** // ** TRACKING ITERATIONS INITIALIZATION ** @@ -895,6 +944,20 @@ InitStatus CbmL1::Init() LOG(info) << " ToF: " << fNTofStations; LOG(info) << " Total: " << fNStations; + // monitor the material + { + LOG(info) << "\033[31;1m-------------------- L1 material -----------------------------\033[0m"; + fMaterialMonitor.clear(); + for (int i = 0; i < fNStations; i++) { + L1MaterialMonitor m(&(fpAlgo->GetParameters()->GetThicknessMaps()[i]), Form("station %d", i)); + LOG(info) << m.ToString(); + fMaterialMonitor.push_back(m); + } + LOG(info) << "\033[31;1m-------------------------------------------------------------\033[0m"; + } + + DumpMaterialToFile("L1material.root"); + return kSUCCESS; } @@ -949,6 +1012,18 @@ void CbmL1::Reconstruct(CbmEvent* event) ReadEvent(event); + // Material monitoring: mark active areas + { + for (L1HitIndex_t i = 0; i < fpAlgo->GetInputData().GetNhits(); i++) { + const L1Hit& h = fpAlgo->GetInputData().GetHit(i); + const L1Station& sta = fpAlgo->GetParameters()->GetStation(h.iSt); + fscal x, y; + std::tie(x, y) = sta.ConvUVtoXY<fscal>(h.u, h.v); + fMaterialMonitor[h.iSt].MarkActiveBin(x, y); + } + } + // + if (fPerformance) { HitMatch(); // calculate the max number of Hits\mcPoints on continuous(consecutive) stations @@ -1051,6 +1126,18 @@ void CbmL1::Reconstruct(CbmEvent* event) // ----- Finish CbmStsFitPerformanceTask task ----------------------------- void CbmL1::Finish() { + + // monitor the material + { + LOG(info) << "\033[31;1m-------------------- L1 material -----------------------------\033[0m"; + + for (int i = 0; i < fNStations; i++) { + LOG(info) << fMaterialMonitor[i].ToString(); + } + + LOG(info) << "\033[31;1m--------------------------------------------------------------\033[0m"; + } + TDirectory* curr = gDirectory; TFile* currentFile = gFile; @@ -1351,50 +1438,6 @@ void CbmL1::WriteSIMDKFData() } } -// --------------------------------------------------------------------------------------------------------------------- -// NOTE: this function should be called before fInitManager.SendParameters(fpAlgo) -std::vector<L1Material> CbmL1::ReadMaterialBudget(L1DetectorID detectorID) -{ - std::vector<L1Material> result {}; - if (fMatBudgetFileName.find(detectorID) != fMatBudgetFileName.end()) { - TFile* oldFile = gFile; - TDirectory* oldDir = gDirectory; - - auto rlFile = TFile(fMatBudgetFileName.at(detectorID)); - if (rlFile.IsZombie()) { LOG(fatal) << "File " << fMatBudgetFileName.at(detectorID) << " is zombie!"; } - else { - LOG(info) << "Reading material budget for " << GetDetectorName(detectorID) << " from file " - << fMatBudgetFileName.at(detectorID); - } - - result.resize(fInitManager.GetNstationsGeometry(detectorID)); - TString stationNamePrefix = "Radiation Thickness [%], Station"; - - // NOTE: Loop over geometry stations. We probably do not know which stations are active/inactive (S.Zharko) - for (int iSt = 0; iSt < fInitManager.GetNstationsGeometry(detectorID); ++iSt) { - // TODO: Unify material table names (S.Zharko) - TString stationName = stationNamePrefix + (detectorID == L1DetectorID::kMvd ? iSt : iSt + 1); - auto* hStaRadLen = rlFile.Get<TProfile2D>(stationName); - LOG_IF(fatal, !hStaRadLen) << "CbmL1: material budget profile " << stationName << " does not exist in file " - << fMatBudgetFileName.at(detectorID); - int nBins = hStaRadLen->GetNbinsX(); - float rMax = hStaRadLen->GetXaxis()->GetXmax(); - result[iSt].SetBins(nBins, rMax); - for (int iBinX = 0; iBinX < nBins; ++iBinX) { - for (int iBinY = 0; iBinY < nBins; ++iBinY) { - result[iSt].SetRadThickBin(iBinX, iBinY, 0.01 * hStaRadLen->GetBinContent(iBinX, iBinY)); - } // iBinX - } // iBinY - LOG(info) << "- station " << iSt; - } // iSt - gFile = oldFile; - gDirectory = oldDir; - } // if mat budget file found - else { - LOG(fatal) << "No material budget file is found for " << GetDetectorName(detectorID); - } - return result; -} // --------------------------------------------------------------------------------------------------------------------- // @@ -1407,3 +1450,36 @@ double CbmL1::boundedGaus(double sigma) } while (fabs(x) > 3.5 * sigma); return x; } + +void CbmL1::DumpMaterialToFile(TString fileName) +{ + TFile* f = new TFile(fileName, "RECREATE"); + f->cd(); + const L1Parameters& param = *fpAlgo->GetParameters(); + const L1MaterialContainer_t& map = param.GetThicknessMaps(); + for (int ist = 0; ist < param.GetNstationsActive(); ist++) { + //const L1Station& st = param.GetStation(ist); + const L1Material& mat = map[ist]; + + TString name = Form("tracking_station%d", ist); + TString title = Form("Tracking station %d: Rad. thickness in [%s]. Z region [%.2f,%.2f] cm.", ist, "%", + mat.GetZmin(), mat.GetZmax()); + + if (fMatBudgetParallelProjection) { title += " Horisontal projection."; } + else { + title += " Radial projection."; + } + + TH2F* h = new TH2F(name.Data(), title.Data(), mat.GetNbins(), -mat.GetXYmax(), mat.GetXYmax(), mat.GetNbins(), + -mat.GetXYmax(), mat.GetXYmax()); + h->GetXaxis()->SetTitle("X [cm]"); + h->GetYaxis()->SetTitle("Y [cm]"); + + for (int ix = 0; ix < mat.GetNbins(); ix++) { + for (int iy = 0; iy < mat.GetNbins(); iy++) { + h->SetBinContent(ix + 1, iy + 1, 100. * mat.GetRadThickBin(ix, iy)); + } + } + } + f->Write(); +} diff --git a/reco/L1/CbmL1.h b/reco/L1/CbmL1.h index afd5bd5dec..492c765cff 100644 --- a/reco/L1/CbmL1.h +++ b/reco/L1/CbmL1.h @@ -59,6 +59,7 @@ #include "L1EventEfficiencies.h" #include "L1IODataManager.h" #include "L1InitManager.h" +#include "L1MaterialMonitor.h" class L1Algo; class L1Event; @@ -226,6 +227,19 @@ public: if (fileName != "") fMatBudgetFileName[L1DetectorID::kTof] = fileName; } + /// Sets material budget binning + void SetMaterialBudgetNbins(int nBinsPerDimension) { fMatBudgetNbins = nBinsPerDimension; } + + /// Sets material budget n rays per dimansion in each bin + void SetMaterialBudgetNrays(int nRaysPerDimension) { fMatBudgetNrays = nRaysPerDimension; } + + /// Sets material budget minimal bin size in cm + void SetMaterialBudgetPitch(double pitchCm) { fMatBudgetPitch = pitchCm; } + + /// Calculate material budget with rays, parallel to Z axis + /// Only needed in debug mode to produce detailed picture of the material + void SetMaterialBudgetParallelProjection() { fMatBudgetParallelProjection = true; } + /// \brief Sets a name for the input configuration file /// If the file is undefined, default tracking parameters will be used. Otherwise, the default parameters will be /// overridden with ones from the configuration file @@ -268,11 +282,6 @@ public: return ""; } - /// Reads material budget information: station thickness in units of radiation length vs. point at the XY plane - /// \param detectorID ID of a detector subsystem - std::vector<L1Material> ReadMaterialBudget(L1DetectorID detectorID); - - void SetExtrapolateToTheEndOfSTS(bool b) { fExtrapolateToTheEndOfSTS = b; } void SetMissingHits(bool value) { fMissingHits = value; } void SetStsOnlyMode() { fTrackingMode = L1Algo::TrackingMode::kSts; } @@ -436,6 +445,8 @@ private: static std::istream& eatwhite(std::istream& is); // skip spaces static void writedir2current(TObject* obj); // help procedure + void DumpMaterialToFile(TString fileName); + private: std::string fInputDataFilename = ""; ///< File name to read/write input hits std::string fsInputSearchWindowsFilename = ""; ///< File name to read search windows @@ -640,11 +651,20 @@ private: std::unordered_map<L1DetectorID, TString> fMatBudgetFileName {}; ///< Map for material budget file names vs. detectorID + int fMatBudgetNbins {100}; ///< n bins in mat budget maps (fMatBudgetNbins x fMatBudgetNbins) + int fMatBudgetNrays {3}; ///< material budget n rays per dimansion in each bin + double fMatBudgetPitch {0.1}; ///< material budget minimal bin size in cm + + bool fMatBudgetParallelProjection {false}; ///< Calculate material budget with rays, parallel to Z axis + ///< Only needed in debug mode to produce detailed picture of the material + bool fExtrapolateToTheEndOfSTS {false}; KFTopoPerformance* fTopoPerformance {nullptr}; L1EventEfficiencies fEventEfficiency {}; // average efficiencies + std::vector<L1MaterialMonitor> fMaterialMonitor {}; ///< material monitoring + ClassDef(CbmL1, 0); }; diff --git a/reco/L1/L1Algo/L1BaseStationInfo.cxx b/reco/L1/L1Algo/L1BaseStationInfo.cxx index 2ac0fa4b52..ae211888c0 100644 --- a/reco/L1/L1Algo/L1BaseStationInfo.cxx +++ b/reco/L1/L1Algo/L1BaseStationInfo.cxx @@ -338,7 +338,6 @@ void L1BaseStationInfo::SetMaterialMap(const L1Material& thicknessMap) { if (!fInitController.GetFlag(EInitKey::kThicknessMap)) { fThicknessMap = thicknessMap; - fThicknessMap.Repare(); fInitController.SetFlag(EInitKey::kThicknessMap); } else { @@ -353,7 +352,6 @@ void L1BaseStationInfo::SetMaterialMap(L1Material&& thicknessMap) noexcept { if (!fInitController.GetFlag(EInitKey::kThicknessMap)) { fThicknessMap = std::move(thicknessMap); - fThicknessMap.Repare(); fInitController.SetFlag(EInitKey::kThicknessMap); } else { diff --git a/reco/L1/L1Algo/L1Material.cxx b/reco/L1/L1Algo/L1Material.cxx index ada685d20d..3393490ebb 100644 --- a/reco/L1/L1Algo/L1Material.cxx +++ b/reco/L1/L1Algo/L1Material.cxx @@ -17,39 +17,6 @@ * L1Material class * ********************/ -L1Material::L1Material() {} - -//------------------------------------------------------------------------------------------------------------------------------------ -// -L1Material::~L1Material() noexcept {} - -//------------------------------------------------------------------------------------------------------------------------------------ -// -L1Material::L1Material(const L1Material& other) - : fNbins(other.fNbins) - , fRmax(other.fRmax) - , fFactor(other.fFactor) - , fTable(other.fTable) -{ -} - -//------------------------------------------------------------------------------------------------------------------------------------ -// - -L1Material& L1Material::operator=(const L1Material& other) -{ - if (this != &other) { - fNbins = other.fNbins; - fRmax = other.fRmax; - fFactor = other.fFactor; - fTable.clear(); - fTable.resize(other.fTable.size()); - for (size_t i = 0; i < other.fTable.size(); ++i) { - fTable[i] = other.fTable[i]; - } - } - return *this; -} //------------------------------------------------------------------------------------------------------------------------------------ // @@ -66,14 +33,23 @@ L1Material& L1Material::operator=(L1Material&& other) noexcept return *this; } +//------------------------------------------------------------------------------------------------------------------------------------ +int L1Material::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 L1Material::GetRadThickScal(float x, float y) const { - x = (x < fRmax && x >= -fRmax) ? x : 0; - y = (y < fRmax && y >= -fRmax) ? y : 0; - int i = static_cast<int>((x + fRmax) * fFactor); - int j = static_cast<int>((y + fRmax) * fFactor); + 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]; @@ -89,13 +65,12 @@ fvec L1Material::GetRadThickVec(fvec x, fvec y) const return r; } - //------------------------------------------------------------------------------------------------------------------------------------ // void L1Material::CheckConsistency() const { /* (i) Check, if the object was initialized */ - if (fNbins < 1 || std::isnan(fRmax)) { throw std::logic_error("L1Material: class was not initialized"); } + if (IsNaN()) { throw std::logic_error("L1Material: class was not initialized"); } /* (ii) Check if the thickness values correct (non-negative) */ bool isThicknessOk = true; @@ -113,24 +88,34 @@ void L1Material::CheckConsistency() const //------------------------------------------------------------------------------------------------------------------------------------ // -void L1Material::SetBins(int nBins, float stationSize) +void L1Material::Initialize(int nBins, float xyMax, float zRef, float zMin, float zMax) { fNbins = nBins; - fRmax = stationSize; + fXYmax = xyMax; + fZref = zRef; + fZmin = zMin; + fZmax = zMax; if (fNbins < 1) { std::stringstream aStream; - aStream << "L1Material: object cannot be initialized with non-positive nBins = " << nBins; + aStream << "L1Material: object cannot be initialized with non-positive nBins = " << fNbins; + throw std::logic_error(aStream.str()); + } + + if (fXYmax < 0.) { + std::stringstream aStream; + aStream << "L1Material: object cannot be initialized with non-positive XYmax = " << fXYmax << " [cm]"; throw std::logic_error(aStream.str()); } - if (stationSize < 0) { + if (!((fZmin <= fZref) && (fZref <= zMax))) { std::stringstream aStream; - aStream << "L1Material: object cannot be initialized with non-positive stationStation = " << stationSize << " [cm]"; + aStream << "L1Material: object cannot be initialized with inconsistent Z: min " << fZmin << " ref " << fZref + << " max " << fZmax << " [cm]"; throw std::logic_error(aStream.str()); } - fFactor = 0.5 * fNbins / fRmax; + fFactor = 0.5 * fNbins / fXYmax; fTable.resize(fNbins * fNbins); } @@ -139,83 +124,10 @@ void L1Material::SetBins(int nBins, float stationSize) void L1Material::Swap(L1Material& other) noexcept { std::swap(fNbins, other.fNbins); - std::swap(fRmax, other.fRmax); + 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) } - -void L1Material::Repare() -{ - // correction of the material map: fill empty bins - // we assume that bins with radiation thickness == 0. are lacking statistics - - const float cut = 1.e-8; - const int n = GetNbins(); - std::vector<float> oldMap(n * n, 0.); - - bool repeat = 1; - - while (repeat) { // until there are empty bins - - oldMap.assign(oldMap.size(), 0.); - for (int iy = 0; iy < n; ++iy) { - for (int ix = 0; ix < n; ++ix) { - oldMap[iy * n + ix] = GetRadThickBin(ix, iy); - } - } - - repeat = 0; - for (int iy = 0; iy < n; ++iy) { - for (int ix = 0; ix < n; ++ix) { - if (oldMap[iy * n + ix] >= cut) continue; - - double sum = 0.; - double sumw = 0.; - // look left - for (int i = ix - 1; i >= 0; --i) { - float v = oldMap[iy * n + i]; - if (v >= cut) { - double w = 1. / (ix - i); - sum += w * v; - sumw += w; - break; - } - } - // look right - for (int i = ix + 1; i < n; ++i) { - float v = oldMap[iy * n + i]; - if (v >= cut) { - double w = 1. / (i - ix); - sum += w * v; - sumw += w; - break; - } - } - // look down - for (int i = iy - 1; i >= 0; --i) { - float v = oldMap[i * n + ix]; - if (v >= cut) { - double w = 1. / (iy - i); - sum += w * v; - sumw += w; - break; - } - } - // look up - for (int i = iy + 1; i < n; ++i) { - float v = oldMap[i * n + ix]; - if (v >= cut) { - double w = 1. / (i - iy); - sum += w * v; - sumw += w; - break; - } - } - if (sumw > 1.e-8) { - SetRadThickBin(ix, iy, sum / sumw); - repeat = true; - } - } - } - } -} diff --git a/reco/L1/L1Algo/L1Material.h b/reco/L1/L1Algo/L1Material.h index 3dcf4924a6..2570fe9065 100644 --- a/reco/L1/L1Algo/L1Material.h +++ b/reco/L1/L1Algo/L1Material.h @@ -20,15 +20,13 @@ class L1Material { public: /// Constructor - /// \param nBins Number of rows or columns - /// \param stationSize Size of station in x and y dimensions [cm] - L1Material(); + L1Material() = default; /// Copy constructor - L1Material(const L1Material& other); + L1Material(const L1Material& other) = default; /// Copy assignment operator - L1Material& operator=(const L1Material& other); + L1Material& operator=(const L1Material& other) = default; /// Move constructor L1Material(L1Material&& other) noexcept; @@ -37,16 +35,32 @@ public: L1Material& operator=(L1Material&& other) noexcept; /// Destructor - ~L1Material() noexcept; + ~L1Material() noexcept = default; /// Gets number of bins (rows or columns) of the material table int GetNbins() const { return fNbins; } + /// Gets radius in cm of the material table + float GetXYmax() const { return fXYmax; } + + /// Gets reference Z of the material in cm + float GetZref() const { return fZref; } + + /// Gets minimal Z of the collected material in cm + float GetZmin() const { return fZmin; } + + /// Gets maximal Z of the collected material in cm + float GetZmax() const { return fZmax; } + /// 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]; } + /// 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]; } + /// 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] @@ -59,7 +73,11 @@ public: cbm::algo::ca::fvec GetRadThickVec(cbm::algo::ca::fvec x, cbm::algo::ca::fvec y) const; /// Checks, if the fields are NaN - bool IsNaN() const { return undef::IsUndefined(fNbins) || undef::IsUndefined(fRmax) || undef::IsUndefined(fFactor); } + bool IsNaN() const + { + return undef::IsUndefined(fNbins) || undef::IsUndefined(fXYmax) || undef::IsUndefined(fFactor) + || undef::IsUndefined(fZref) || undef::IsUndefined(fZmin) || undef::IsUndefined(fZmax); + } /// Verifies class invariant consistency void CheckConsistency() const; @@ -75,19 +93,22 @@ public: /// 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 stationSize Size of station in x and y dimensions [cm] - void SetBins(int nBins, float stationSize); + /// \param xyMax Size of station in x and y dimensions [cm] + void Initialize(int nBins, float xyMax, float zRef, float zMin, float zMax); /// Swap method void Swap(L1Material& other) noexcept; - /// repare the map - fill empty bins - void Repare(); + /// Get bin index for (x,y). Returns -1 when outside of the map + int GetBin(float x, float y) const; private: - int fNbins = undef::kI32; ///< Number of rows (columns) in the material budget table - float fRmax = undef::kF; ///< Size of the station in x and y dimensions [cm] + int fNbins = undef::kI32; ///< Number of rows (== N columns) in the material budget table + float fXYmax = undef::kF; ///< Size of the station in x and y dimensions [cm] float fFactor = undef::kF; ///< Factor used in the recalculation of point coordinates to row/column id + float fZref = undef::kF; ///< Reference Z of the collected material [cm] + float fZmin = undef::kF; ///< Minimal Z of the collected material [cm] + float fZmax = undef::kF; ///< Minimal Z of the collected material [cm] std::vector<float> fTable {}; ///< Material budget table /// Serialization function @@ -96,8 +117,11 @@ private: void serialize(Archive& ar, const unsigned int) { ar& fNbins; - ar& fRmax; + ar& fXYmax; ar& fFactor; + ar& fZref; + ar& fZmin; + ar& fZmax; ar& fTable; } } _fvecalignment; diff --git a/reco/L1/L1Algo/L1MaterialMonitor.cxx b/reco/L1/L1Algo/L1MaterialMonitor.cxx new file mode 100644 index 0000000000..f42567447a --- /dev/null +++ b/reco/L1/L1Algo/L1MaterialMonitor.cxx @@ -0,0 +1,146 @@ +/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt + SPDX-License-Identifier: GPL-3.0-only + Authors: Sergey Gorbunov [committer] */ + +#include "L1MaterialMonitor.h" + +#include <Logger.h> + +#include <iomanip> +#include <sstream> + + +void L1MaterialMonitor::SetMaterial(const L1Material* materialMap) +{ + /// construction + /// \param materialMap Matareial map to be monitored + + fMaterial = materialMap; + + fActiveBinMap.resize(0); + + if (fMaterial) { fActiveBinMap.assign(fMaterial->GetNbins() * fMaterial->GetNbins(), 0); } + + fNhitsTotal = 0; + fNhitsOutside = 0; + + EvaluateStatistics(); +} + +void L1MaterialMonitor::MarkActiveBin(float x, float y) +{ + /// mark a bin as active + if (!fMaterial) { + LOG(fatal) << "L1MaterialMonitor: material map is not set"; + return; + } + int i = fMaterial->GetBin(x, y); + fNhitsTotal++; + if (i < 0) { fNhitsOutside++; } + else { + fActiveBinMap[i] = 1; + } +} + +void L1MaterialMonitor::EvaluateStatistics() +{ + /// update values of statistical variables with respect to the active map + + fActiveNbins = 0; + fPassiveNbins = 0; + + fActiveRadThickMin = undef::kD; + fActiveRadThickMax = undef::kD; + fActiveRadThickMean = undef::kD; + + fPassiveRadThickMin = undef::kD; + fPassiveRadThickMax = undef::kD; + fPassiveRadThickMean = undef::kD; + + if (!fMaterial) { return; } + + int nBins = fMaterial->GetNbins() * fMaterial->GetNbins(); + + if (nBins != (int) fActiveBinMap.size()) { + LOG(fatal) << "L1MaterialMonitor: map of active bins is not consistent with the mterial map: nbins " + << fActiveBinMap.size() << " != " << nBins; + return; + } + + for (int i = 0; i < nBins; i++) { + double r = fMaterial->GetRadThickBin(i); + if (fActiveBinMap[i]) { // active material + if (fActiveNbins == 0) { + fActiveRadThickMin = r; + fActiveRadThickMax = r; + fActiveRadThickMean = r; + } + else { + fActiveRadThickMin = std::min(fActiveRadThickMin, r); + fActiveRadThickMax = std::max(fActiveRadThickMax, r); + fActiveRadThickMean += r; + } + fActiveNbins++; + } + else { + // passive material + if (fPassiveNbins == 0) { + fPassiveRadThickMin = r; + fPassiveRadThickMax = r; + fPassiveRadThickMean = r; + } + else { + fPassiveRadThickMin = std::min(fPassiveRadThickMin, r); + fPassiveRadThickMax = std::max(fPassiveRadThickMax, r); + fPassiveRadThickMean += r; + } + fPassiveNbins++; + } + } + if (fActiveNbins + fPassiveNbins != nBins) { + LOG(fatal) << "L1MaterialMonitor: wrong calculation of N passive / active bins "; + } + if (fActiveNbins > 0) { fActiveRadThickMean /= fActiveNbins; } + if (fPassiveNbins > 0) { fPassiveRadThickMean /= fPassiveNbins; } +} + + +TString L1MaterialMonitor::ToString() +{ + /// print statistics to a string + + EvaluateStatistics(); + + TString s(Form("material map %s. ", fName.Data())); + + if (fActiveNbins > 0) { + s += Form("Active material RL: min %.2f%%, max %.2f%%, mean %.2f%%. ", fActiveRadThickMin * 100., + fActiveRadThickMax * 100., fActiveRadThickMean * 100); + } + else { + if (fNhitsTotal > 0) { s += "No active material. "; } + else { + s += "No hits to identify active areas. "; + } + } + + if (fPassiveNbins > 0) { + s += Form("Passive material RL: min %.2f%%, max %.2f%%, mean %.2f%%. ", fPassiveRadThickMin * 100., + fPassiveRadThickMax * 100., fPassiveRadThickMean * 100); + } + else { + s += "No passive material. "; + } + + if (fNhitsTotal > 0) { + if (fNhitsOutside == 0) { s += "No hits outside of the map. "; } + else { + s += Form("There are %f%% hits outside the map!!! ", (100. * fNhitsOutside) / fNhitsTotal); + } + } + else { + s += "No hit statistics. "; + } + + return s; +} diff --git a/reco/L1/L1Algo/L1MaterialMonitor.h b/reco/L1/L1Algo/L1MaterialMonitor.h new file mode 100644 index 0000000000..63cf63c059 --- /dev/null +++ b/reco/L1/L1Algo/L1MaterialMonitor.h @@ -0,0 +1,105 @@ +/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt + SPDX-License-Identifier: GPL-3.0-only + Authors: Sergey Gorbunov [committer] */ + +#ifndef L1MaterialMonitor_h +#define L1MaterialMonitor_h + + +#include "TString.h" + +#include <vector> + +#include "L1Material.h" + + +/// Class to collect statistics for L1Material +/// +class L1MaterialMonitor { +public: + /// default constructor + L1MaterialMonitor() : L1MaterialMonitor(nullptr) {} + + /// constructor + /// \param materialMap Matareial map to be monitored + /// \param name Name of the material map + L1MaterialMonitor(const L1Material* materialMap, TString name = "") : fName(name) { SetMaterial(materialMap); } + + /// construction + /// \param materialMap Matareial map to be monitored + void SetMaterial(const L1Material* materialMap); + + /// construction + /// \param name Name of the material map + void SetName(TString name) { fName = name; } + + /// reset the map of active bins + void ResetActiveBins() + { + for (auto& v : fActiveBinMap) { + v = 0; + } + } + + /// mark a bin as active + void MarkActiveBin(float x, float y); + + /// update values of statistical variables with respect to the active map + void EvaluateStatistics(); + + /// print statistics to a string + TString ToString(); + + /// get number of active bins in the map + int GetActiveNbins() const { return fActiveNbins; } + + /// get minimal radiation thickness among all active bins + double GetActiveRadThickMin() const { return fPassiveRadThickMin; } + + /// get maximal radiation thickness among all active bins + double GetActiveRadThickMax() const { return fPassiveRadThickMax; } + + /// get average radiation thickness among all active bins + double GetActiveRadThickMean() const { return fPassiveRadThickMean; } + + + /// get number of passive bins in the map + int GetPassiveNbins() const { return fPassiveNbins; } + + /// get minimal radiation thickness among all passive bins + double GetPassiveRadThickMin() const { return fPassiveRadThickMin; } + + /// get maximal radiation thickness among all passive bins + double GetPassiveRadThickMax() const { return fPassiveRadThickMax; } + + /// 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); } + + /// get the number of processed hits + double GetNhits() const { return fNhitsTotal; } + +private: + const L1Material* fMaterial {nullptr}; ///< Pointer to L1Material + TString fName {}; ///< Name of the material map + + std::vector<char> fActiveBinMap {}; ///< Map of active bins in the material map (bins where hits appear) + + int fActiveNbins {undef::kI32}; ///< Active material: number of bins + double fActiveRadThickMin {undef::kD}; ///< Active material: minimal thickness + double fActiveRadThickMax {undef::kD}; ///< Active material: maximal thickness + double fActiveRadThickMean {undef::kD}; ///< Active material: average thickness + + int fPassiveNbins {undef::kI32}; ///< Passive material: number of bins + double fPassiveRadThickMin {undef::kD}; ///< Passive material: minimal thickness + double fPassiveRadThickMax {undef::kD}; ///< Passive material: maximal thickness + double fPassiveRadThickMean {undef::kD}; ///< Passive material: average thickness + + unsigned long fNhitsTotal {0}; ///< number of hits in statistics + unsigned long fNhitsOutside {0}; ///< number of hits outside the material map +}; + +#endif diff --git a/reco/L1/L1LinkDef.h b/reco/L1/L1LinkDef.h index 5054f2b2b1..7ac1d56817 100644 --- a/reco/L1/L1LinkDef.h +++ b/reco/L1/L1LinkDef.h @@ -44,5 +44,6 @@ //#pragma link C++ class cbm::ca::IdealHitProducer < L1DetectorID::kTrd > + ; //#pragma link C++ class cbm::ca::IdealHitProducer < L1DetectorID::kTof > + ; #pragma link C++ class cbm::ca::IdealHitProducer + ; +#pragma link C++ class ca::tools::MaterialHelper + ; #endif diff --git a/reco/L1/catools/CaToolsMaterialHelper.cxx b/reco/L1/catools/CaToolsMaterialHelper.cxx new file mode 100644 index 0000000000..e2057b096c --- /dev/null +++ b/reco/L1/catools/CaToolsMaterialHelper.cxx @@ -0,0 +1,304 @@ +/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt + SPDX-License-Identifier: GPL-3.0-only + Authors: Sergey Gorbunov [committer] */ + + +#include "CaToolsMaterialHelper.h" + +#include "Logger.h" + +#include "TGeoManager.h" +#include "TGeoMedium.h" +#include "TGeoNavigator.h" +#include "TGeoNode.h" + +#include <iostream> +#include <thread> +using namespace ca::tools; + +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) << "------------------------------------------------------------------------------------"; + + InitThreads(); +} + + +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(); +} + + +L1Material 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; + } + + L1Material 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) { + std::cout << "ray at " << x << " " << y << " " << z << " zMin " << zMin << " zMax " << zMax << std::endl; + } + + TGeoNode* node = navi->InitTrack(x, y, z, tx / t, ty / t, 1. / t); + + bool doContinue = 1; + for (int iStep = 0; doContinue; iStep++) { + + if (!node) { LOG(fatal) << "ca::tools::MaterialHelper: TGeoNavigator can not find the geo node"; } + + 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) { + std::cout << " RL " << radLen << " next pt " << navi->GetCurrentPoint()[0] << " " + << navi->GetCurrentPoint()[1] << " " << navi->GetCurrentPoint()[2] << std::endl; + } + + double zNew = navi->GetCurrentPoint()[2]; + if (zNew >= zMax) { + zNew = zMax; + doContinue = 0; + } + + //if (zNew < z) { + //std::cout << " z curr " << z << " z new " << zNew << " dz " << zNew - z << std::endl; } + 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) { std::cout << " radThick " << radThick << std::endl; } + //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 std::move(matBudget); +} diff --git a/reco/L1/catools/CaToolsMaterialHelper.h b/reco/L1/catools/CaToolsMaterialHelper.h new file mode 100644 index 0000000000..487df8679a --- /dev/null +++ b/reco/L1/catools/CaToolsMaterialHelper.h @@ -0,0 +1,95 @@ +/* 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 "Rtypes.h" +#include "TObject.h" + +#include "L1Material.h" + +class TGeoNavigator; + +namespace ca +{ + namespace tools + { + /// Class ca::tools::Debugger helps to debug CA tracking + /// + /// + /// class to create L1Material 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 + /// + L1Material GenerateMaterialMap(double zRef, double zMin, double zMax, double xyMax, int nBinsDim); + + 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 + + ClassDef(MaterialHelper, 0); + }; + + } // namespace tools +} // namespace ca + + +#endif -- GitLab