From 5d743208438375a4eff1a9fd14b7b262435c84e1 Mon Sep 17 00:00:00 2001 From: "s.zharko@gsi.de" <s.zharko@gsi.de> Date: Thu, 9 May 2024 00:17:41 +0200 Subject: [PATCH] 1) the input and output online CA QA are merged, the amount of histograms was reduced 2) the mcbm_beam_2024_05_08_nickel setup is added to the CbmMcbmUtils (for RunID > 2917) --- algo/CMakeLists.txt | 8 +- algo/ca/TrackingChain.cxx | 43 +-- algo/ca/TrackingChain.h | 7 +- algo/ca/core/utils/CaTrackingMonitor.h | 7 +- algo/ca/qa/CaInputQa.cxx | 272 ----------------- algo/ca/qa/CaInputQa.h | 91 ------ algo/ca/qa/CaOutputQa.cxx | 178 ----------- algo/ca/qa/CaOutputQa.h | 63 ---- algo/ca/qa/CaQa.cxx | 402 +++++++++++++++++++++++++ algo/ca/qa/CaQa.h | 156 ++++++++++ algo/ca/qa/CaQaBase.cxx | 36 --- algo/ca/qa/CaQaBase.h | 100 ------ algo/qa/PadConfig.h | 16 + core/base/utils/CbmMcbmUtils.cxx | 3 +- 14 files changed, 594 insertions(+), 788 deletions(-) delete mode 100644 algo/ca/qa/CaInputQa.cxx delete mode 100644 algo/ca/qa/CaInputQa.h delete mode 100644 algo/ca/qa/CaOutputQa.cxx delete mode 100644 algo/ca/qa/CaOutputQa.h create mode 100644 algo/ca/qa/CaQa.cxx create mode 100644 algo/ca/qa/CaQa.h delete mode 100644 algo/ca/qa/CaQaBase.cxx delete mode 100644 algo/ca/qa/CaQaBase.h diff --git a/algo/CMakeLists.txt b/algo/CMakeLists.txt index bc7874f883..b84c480fa2 100644 --- a/algo/CMakeLists.txt +++ b/algo/CMakeLists.txt @@ -157,9 +157,7 @@ set(SRCS qa/unpack/StsDigiQa.cxx ca/TrackingSetup.cxx ca/TrackingChain.cxx - ca/qa/CaInputQa.cxx - ca/qa/CaOutputQa.cxx - ca/qa/CaQaBase.cxx + ca/qa/CaQa.cxx ) set(BUILD_INFO_CXX ${CMAKE_CURRENT_BINARY_DIR}/base/BuildInfo.cxx) @@ -284,9 +282,7 @@ install( ca/TrackingChainConfig.h # NOTE: SZh 20.11.2023: # The ca/qa directory depends on the online qa classes, so for now it has to be a part of the Algo library. - ca/qa/CaInputQa.h - ca/qa/CaOutputQa.h - ca/qa/CaQaBase.h + ca/qa/CaQa.h DESTINATION include/ ) diff --git a/algo/ca/TrackingChain.cxx b/algo/ca/TrackingChain.cxx index ee2a90edf7..39178f8573 100644 --- a/algo/ca/TrackingChain.cxx +++ b/algo/ca/TrackingChain.cxx @@ -33,20 +33,15 @@ using cbm::algo::ca::EDetectorID; using cbm::algo::ca::Framework; using cbm::algo::ca::HitTypes_t; using cbm::algo::ca::InitManager; -using cbm::algo::ca::InputQa; -using cbm::algo::ca::OutputQa; using cbm::algo::ca::Parameters; +using cbm::algo::ca::Qa; using cbm::algo::ca::Track; using cbm::algo::ca::constants::clrs::CL; // clear text using cbm::algo::ca::constants::clrs::GNb; // grin bald text // --------------------------------------------------------------------------------------------------------------------- // -TrackingChain::TrackingChain(std::shared_ptr<HistogramSender> histoSender) - : fInputQa(InputQa(histoSender)) - , fOutputQa(OutputQa(histoSender)) -{ -} +TrackingChain::TrackingChain(std::shared_ptr<HistogramSender> histoSender) : fQa(Qa(histoSender)) {} // --------------------------------------------------------------------------------------------------------------------- // @@ -94,14 +89,9 @@ void TrackingChain::Init() fCaFramework.Init(ca::Framework::TrackingMode::kMcbm); // ------ Initialize QA modules - if (fInputQa.IsSenderDefined()) { - fInputQa.RegisterParameters(&fCaFramework.GetParameters()); - fInputQa.Init(); - } - - if (fOutputQa.IsSenderDefined()) { - fOutputQa.RegisterParameters(&fCaFramework.GetParameters()); - fOutputQa.Init(); + if (fQa.IsSenderDefined()) { + fQa.RegisterParameters(&fCaFramework.GetParameters()); + fQa.Init(); } } @@ -208,22 +198,13 @@ TrackingChain::Output_t TrackingChain::PrepareOutput() // QA - if (fInputQa.IsSenderDefined()) { - fCaMonitorData.StartTimer(ca::ETimer::InputQa); - fInputQa.RegisterInputData(&fCaFramework.GetInputData()); - fInputQa.RegisterTracks(&output.tracks); - fInputQa.RegisterRecoHitIndices(&fCaFramework.fRecoHits); - fInputQa.Exec(); - fCaMonitorData.StopTimer(ca::ETimer::InputQa); - } - - if (fOutputQa.IsSenderDefined()) { - fCaMonitorData.StartTimer(ca::ETimer::OutputQa); - fOutputQa.RegisterInputData(&fCaFramework.GetInputData()); - fOutputQa.RegisterTracks(&output.tracks); - fOutputQa.RegisterRecoHitIndices(&fCaFramework.fRecoHits); - fOutputQa.Exec(); - fCaMonitorData.StopTimer(ca::ETimer::OutputQa); + if (fQa.IsSenderDefined()) { + fCaMonitorData.StartTimer(ca::ETimer::Qa); + fQa.RegisterInputData(&fCaFramework.GetInputData()); + fQa.RegisterTracks(&output.tracks); + fQa.RegisterRecoHitIndices(&fCaFramework.fRecoHits); + fQa.Exec(); + fCaMonitorData.StopTimer(ca::ETimer::Qa); } fCaMonitorData.StopTimer(ca::ETimer::TrackingChain); diff --git a/algo/ca/TrackingChain.h b/algo/ca/TrackingChain.h index d9e852a732..e96735b5be 100644 --- a/algo/ca/TrackingChain.h +++ b/algo/ca/TrackingChain.h @@ -11,8 +11,7 @@ #include "CaDataManager.h" #include "CaFramework.h" -#include "CaInputQa.h" -#include "CaOutputQa.h" +#include "CaQa.h" #include "CaTrack.h" #include "CaTrackingMonitor.h" #include "CaVector.h" @@ -112,9 +111,7 @@ namespace cbm::algo ca::TrackingMonitorData fCaMonitorData{}; ///< CA monitor data object ca::Framework fCaFramework{}; ///< CA framework instance ca::DataManager fCaDataManager{}; ///< CA data manager - - ca::InputQa fInputQa; ///< CA input QA builder - ca::OutputQa fOutputQa; ///< CA output QA builder + ca::Qa fQa; ///< CA QA builder std::shared_ptr<TrackingSetup> fpSetup = nullptr; ///< setup interface diff --git a/algo/ca/core/utils/CaTrackingMonitor.h b/algo/ca/core/utils/CaTrackingMonitor.h index 0ba75e3459..da37a1632a 100644 --- a/algo/ca/core/utils/CaTrackingMonitor.h +++ b/algo/ca/core/utils/CaTrackingMonitor.h @@ -70,8 +70,7 @@ namespace cbm::algo::ca MergeClones, StoreTracksWindow, StoreTracksFinal, - InputQa, - OutputQa, + Qa, kEND }; @@ -103,6 +102,7 @@ namespace cbm::algo::ca SetCounterName(ECounter::UndefinedTrdHit, "undefined TRD hits"); SetCounterName(ECounter::UndefinedTofHit, "undefined TOF hits"); + SetTimerName(ETimer::TrackingChain, "tracking chain"); SetTimerName(ETimer::PrepareInputData, "input data preparation"); SetTimerName(ETimer::Tracking, "algorithm execution"); SetTimerName(ETimer::PrepareTimeslice, "timeslice preparation"); @@ -122,8 +122,7 @@ namespace cbm::algo::ca SetTimerName(ETimer::MergeClones, "clone merger"); SetTimerName(ETimer::StoreTracksWindow, "track storing in window"); SetTimerName(ETimer::StoreTracksFinal, "final track storing"); - SetTimerName(ETimer::InputQa, "input QA"); - SetTimerName(ETimer::OutputQa, "output QA"); + SetTimerName(ETimer::Qa, "QA"); SetRatioKeys({ECounter::TrackingCall, ECounter::SubTS, ECounter::RecoTrack}); } diff --git a/algo/ca/qa/CaInputQa.cxx b/algo/ca/qa/CaInputQa.cxx deleted file mode 100644 index b0c40e0489..0000000000 --- a/algo/ca/qa/CaInputQa.cxx +++ /dev/null @@ -1,272 +0,0 @@ -/* Copyright (C) 2023-2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt - SPDX-License-Identifier: GPL-3.0-only - Authors: Sergei Zharko [committer] */ - -/// \file CaInputQa.cxx -/// \date 01.03.2024 -/// \brief An input QA module for CA tracking (source) -/// \author Sergei Zharko <s.zharko@gsi.de> - -#include "CaInputQa.h" - -#include "CaConstants.h" -#include "CaInputData.h" -#include "CaParameters.h" -#include "TrackingDefs.h" - -#include <algorithm> -#include <fstream> -#include <limits> - -#include <fmt/format.h> - -using cbm::algo::ca::Hit; -using cbm::algo::ca::InputQa; -using cbm::algo::ca::constants::math::Pi; - - -// --------------------------------------------------------------------------------------------------------------------- -// -void InputQa::Init() -{ - using cbm::algo::qa::CanvasConfig; - using cbm::algo::qa::Data; - using cbm::algo::qa::H1D; - using cbm::algo::qa::H2D; - using cbm::algo::qa::PadConfig; - using cbm::algo::qa::Prof2D; - using fmt::format; - - if (!fpSender.get()) { - return; - } - - if (!fpParameters) { - L_(error) << "ca::InputQa::Init(): parameters object is not initialized"; - assert(false); - } - - int nSt = fpParameters->GetNstationsActive(); // Number of active tracking stations - - // ------ Init histograms - - // Occupancy distributions - fvphHitOccupXY.resize(nSt); - fvphHitOccupZX.resize(nSt); - fvphHitOccupZY.resize(nSt); - fvphHitUsageXY.resize(nSt); - - // Station sizes in transverse plane - std::vector<double> vMinX(nSt); - std::vector<double> vMaxX(nSt); - std::vector<double> vMinY(nSt); - std::vector<double> vMaxY(nSt); - - int nBinsXY = 200; - int nBinsZ = 400; - - for (int iSt = 0; iSt < nSt; ++iSt) { - const auto& station = fpParameters->GetStation(iSt); - vMaxX[iSt] = station.GetXmax<double>(); - vMinX[iSt] = station.GetXmin<double>(); - vMaxY[iSt] = station.GetYmax<double>(); - vMinY[iSt] = station.GetYmin<double>(); - double dy = (vMaxY[iSt] - vMinY[iSt]) * kXYZMargin; - double dx = (vMaxX[iSt] - vMinX[iSt]) * kXYZMargin; - vMaxX[iSt] += dx; - vMinX[iSt] -= dx; - vMaxY[iSt] += dy; - vMinY[iSt] -= dy; - } - // Station max - double xMinA = *std::min_element(vMinX.begin(), vMinX.end()); - double xMaxA = *std::max_element(vMaxX.begin(), vMaxX.end()); - double yMinA = *std::min_element(vMinY.begin(), vMinY.end()); - double yMaxA = *std::max_element(vMaxY.begin(), vMaxY.end()); - double zMinA = fpParameters->GetStation(0).GetZ<double>(); - double zMaxA = fpParameters->GetStation(nSt - 1).GetZ<double>(); - { - double dz = (zMaxA - zMinA) * kXYZMargin; - zMinA -= dz; - zMaxA += dz; - } - for (int iSt = 0; iSt < nSt; ++iSt) { - int iStGeo = fpParameters->GetActiveToGeoMap()[iSt]; - auto [detID, iStLoc] = fpParameters->GetGeoToLocalIdMap()[iStGeo]; - for (auto hitSet : kHitSets) { - auto setNm = EHitSet::Input == hitSet ? "input" : "used"; - auto setTl = EHitSet::Input == hitSet ? "Input" : "Used"; - { - auto name = format("hit_{}_occup_xy_sta_{}", setNm, iSt); - auto titl = format("{} hit occupancy in XY plane for station {} ({}{});x [cm];y [cm]", setTl, iSt, - kDetName[detID], iStLoc); - fvphHitOccupXY[iSt][hitSet] = - fQaData.MakeObj<H2D>(name, titl, nBinsXY, vMinX[iSt], vMaxX[iSt], nBinsXY, vMinY[iSt], vMaxY[iSt]); - } - { - auto name = format("hit_{}_occup_zx_sta_{}", setNm, iSt); - auto titl = format("{} hit occupancy in ZX plane for station {} ({}{});z [cm];x [cm]", setTl, iSt, - kDetName[detID], iStLoc); - fvphHitOccupZX[iSt][hitSet] = fQaData.MakeObj<H2D>(name, titl, nBinsZ, zMinA, zMaxA, nBinsXY, xMinA, xMaxA); - } - { - auto name = format("hit_{}_occup_zy_sta_{}", setNm, iSt); - auto titl = format("{} hit occupancy in ZY plane for station {} ({}{});z [cm];y [cm]", setTl, iSt, - kDetName[detID], iStLoc); - fvphHitOccupZY[iSt][hitSet] = fQaData.MakeObj<H2D>(name, titl, nBinsZ, zMinA, zMaxA, nBinsXY, yMinA, yMaxA); - } - } - { - auto name = format("hit_usage_xy_sta_{}", iSt); - auto titl = format("Hit usage in XY plane for station {} ({}{});x [cm];y [cm]", iSt, kDetName[detID], iStLoc); - fvphHitUsageXY[iSt] = - fQaData.MakeObj<Prof2D>(name, titl, nBinsXY, vMinX[iSt], vMaxX[iSt], nBinsXY, vMinY[iSt], vMaxY[iSt], 0., 1.); - } - } - for (auto hitSet : kHitSets) { - constexpr int NBins = 1000000; - auto setNm = EHitSet::Input == hitSet ? "input" : "used"; - auto setTl = EHitSet::Input == hitSet ? "Input" : "Used"; - { - auto name = format("hit_{}_front_key_index", setNm); - auto titl = format("{} hit front key index;ID_{{key}}/N_{{keys}};Count", setTl); - fvphHitFrontKeyIndex[hitSet] = fQaData.MakeObj<H1D>(name, titl, NBins, 0., 1.); - } - { - auto name = format("hit_{}_back_key_index", setNm); - auto titl = format("{} hit back key index;ID_{{key}}/N_{{keys}};Count", setTl); - fvphHitBackKeyIndex[hitSet] = fQaData.MakeObj<H1D>(name, titl, NBins, 0., 1.); - } - } - - // Hit time distributions - fvphHitTime.resize(nSt + 1); // [nSt] - total over stations - for (int iSt = 0; iSt < nSt + 1; ++iSt) { - auto staNm = iSt == nSt ? "" : format("_sta_{}", iSt); - auto staTl = iSt == nSt ? "" : format(" in station {}", iSt); - for (auto hitSet : kHitSets) { - auto setNm = EHitSet::Input == hitSet ? "input" : "used"; - auto setTl = EHitSet::Input == hitSet ? "Input" : "Used"; - auto name = format("hit_{}_rel_time{}", setNm, staNm); - auto titl = format("{} hit relative time{}; #delta t_{{hit}};Count", setTl, staTl); - fvphHitTime[iSt][hitSet] = fQaData.MakeObj<H1D>(name, titl, 10000, -0.1, 1.1); - } - } - - // ----- Init canvases - // Hit occupancies - for (auto hitSet : kHitSets) { - auto setNm = EHitSet::Input == hitSet ? "input" : "used"; - auto setTl = EHitSet::Input == hitSet ? "Input" : "Used"; - { // XY - auto name = format("ca_hit_{}_occupancy_xy", setNm); - auto titl = format("{} hit occupancy in different stations in XY plane", setNm); - auto canv = CanvasConfig(name, titl); - for (int iSt = 0; iSt < nSt; ++iSt) { - auto pad = PadConfig(); - pad.RegisterHistogram(fvphHitOccupXY[iSt][hitSet], "colz"); - canv.AddPadConfig(pad); - } - fQaData.AddCanvasConfig(canv); - } - { // ZX and ZY - auto name = format("ca_hit_{}_occupancy_zx_zy", setNm); - auto titl = format("{} hit occupancy in different stations in ZX and ZY planes", setTl); - auto canv = CanvasConfig(name, titl); - { // ZX - auto pad = PadConfig(); - for (int iSt = 0; iSt < nSt; ++iSt) { - pad.RegisterHistogram(fvphHitOccupZX[iSt][hitSet], (iSt == 0 ? "colz" : "cols same")); - } - canv.AddPadConfig(pad); - } - { // ZY - auto pad = PadConfig(); - for (int iSt = 0; iSt < nSt; ++iSt) { - pad.RegisterHistogram(fvphHitOccupZY[iSt][hitSet], (iSt == 0 ? "colz" : "cols same")); - } - canv.AddPadConfig(pad); - } - fQaData.AddCanvasConfig(canv); - } - } - // Hit usage profiles - { - auto name = format("ca_hit_usage_xy"); - auto titl = format("Hit usage in different stations in XY plane"); - auto canv = CanvasConfig(name, titl); - for (int iSt = 0; iSt < nSt; ++iSt) { - auto pad = PadConfig(); - pad.RegisterHistogram(fvphHitUsageXY[iSt], "colz"); - canv.AddPadConfig(pad); - } - fQaData.AddCanvasConfig(canv); - } - - fQaData.Init(fpSender); -} - -// --------------------------------------------------------------------------------------------------------------------- -// -void InputQa::Exec() -{ - if (!fpSender.get()) { - return; - } - - if (!CheckInit()) { - L_(fatal) << "ca::OutputQa: instance is not initialized"; - assert(false); - } - - int nHitsInput = fpInputData->GetNhits(); - // Map: if hit used in tracking - std::vector<unsigned char> vbHitUsed(nHitsInput, false); - for (int iH : (*fpvRecoHits)) { - vbHitUsed[iH] = true; - } - - // Calculate max and min hit time - const auto& hits = fpInputData->GetHits(); - auto [minTimeIt, maxTimeIt] = - std::minmax_element(hits.cbegin(), hits.cend(), [](const auto& h1, const auto& h2) { return h1.T() < h2.T(); }); - fMinHitTime = minTimeIt->T(); - fMaxHitTime = maxTimeIt->T(); - - // Fill input hit histograms - { - for (int iH = 0; iH < nHitsInput; ++iH) { - const auto& hit = fpInputData->GetHit(iH); - FillHitDistributionsForHitSet(EHitSet::Input, hit); - if (vbHitUsed[iH]) { - FillHitDistributionsForHitSet(EHitSet::Used, hit); - } - int iSt = hit.Station(); - double x = hit.X(); - double y = hit.Y(); - fvphHitUsageXY[iSt]->Fill(x, y, static_cast<double>(vbHitUsed[iH])); - } - } - - fQaData.Send(fpSender); -} - -// --------------------------------------------------------------------------------------------------------------------- -// -void InputQa::FillHitDistributionsForHitSet(InputQa::EHitSet hitSet, const Hit& hit) -{ - int nSt = fpParameters->GetNstationsActive(); - int iSt = hit.Station(); - double x = hit.X(); - double y = hit.Y(); - double z = hit.Z(); - fvphHitOccupXY[iSt][hitSet]->Fill(x, y); - fvphHitOccupZX[iSt][hitSet]->Fill(z, x); - fvphHitOccupZY[iSt][hitSet]->Fill(z, y); - auto nKeys = static_cast<double>(fpInputData->GetNhitKeys()); - fvphHitFrontKeyIndex[hitSet]->Fill(hit.FrontKey() / nKeys); - fvphHitBackKeyIndex[hitSet]->Fill(hit.BackKey() / nKeys); - double relTime = (hit.T() - fMinHitTime) / (fMaxHitTime - fMinHitTime); - fvphHitTime[iSt][hitSet]->Fill(relTime); - fvphHitTime[nSt][hitSet]->Fill(relTime); -} diff --git a/algo/ca/qa/CaInputQa.h b/algo/ca/qa/CaInputQa.h deleted file mode 100644 index adcf82152d..0000000000 --- a/algo/ca/qa/CaInputQa.h +++ /dev/null @@ -1,91 +0,0 @@ -/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt - SPDX-License-Identifier: GPL-3.0-only - Authors: Sergei Zharko [committer] */ - -/// \file CaInputQa.h -/// \date 01.03.2024 -/// \brief A QA module for CA tracking input data (header) -/// \author Sergei Zharko <s.zharko@gsi.de> - -#pragma once - -#include "CaEnumArray.h" -#include "CaQaBase.h" - -namespace cbm::algo::ca -{ - /// \class cbm::algo::ca::qa::InputQa - /// \brief InputQa class for the CA tracking QA (header) - /// - class InputQa : public QaBase { - /// \brief Hit set entries - enum class EHitSet - { - Input, ///< Input hits - Used, ///< Hits used in tracks - kEND - }; - - /// \brief Definition of enum array over EHitSet entries - template<typename T> - using HitSetArray_t = EnumArray<EHitSet, T>; - - /// \brief Array of EHitSet entries for iteration - static constexpr HitSetArray_t<EHitSet> kHitSets = {EHitSet::Input, EHitSet::Used}; - - public: - /// \brief Default destructor - /// \param pSender Pointer to the histogram sender - InputQa(std::shared_ptr<HistogramSender> pSender) : QaBase(pSender, "CA/Input"){}; - - /// \brief Constructor from the configuration object - /// \param config QA configuration object - InputQa() = default; - - /// \brief Copy constructor - InputQa(const InputQa&) = delete; - - /// \brief Move constructor - InputQa(InputQa&&) = delete; - - /// \brief Destructor - ~InputQa() = default; - - /// \brief Copy assignment operator - InputQa& operator=(const InputQa&) = delete; - - /// \brief Move assignment operator - InputQa& operator=(InputQa&&) = delete; - - /// \brief QA execution function - void Exec(); - - /// \brief Initializes the QA - void Init(); - - private: - /// \brief Fills hit distributions - /// \param hitSet Hit set enum entry - /// \param hit Reference to hit - void FillHitDistributionsForHitSet(EHitSet hitSet, const ca::Hit& hit); - - static constexpr double kXYZMargin = 0.05; ///< Margin for occupancy distributions in XY plane - static constexpr int knHitSets = 2; ///< Number of hit sets: input/used - - double fMinHitTime = std::numeric_limits<double>::max(); - double fMaxHitTime = std::numeric_limits<double>::lowest(); - - // Hit distributions vs. station - using OccupHistContainer_t = std::vector<HitSetArray_t<qa::H2D*>>; - OccupHistContainer_t fvphHitOccupXY; ///< hist: Hit occupancy in different stations in XY plane - OccupHistContainer_t fvphHitOccupZX; ///< hist: Hit occupancy in different stations in ZX plane - OccupHistContainer_t fvphHitOccupZY; ///< hist: Hit occupancy in different stations in ZY plane - - std::vector<qa::Prof2D*> fvphHitUsageXY; ///< prof: Hit usage in different stations in XY plane - - HitSetArray_t<qa::H1D*> fvphHitFrontKeyIndex = {nullptr, nullptr}; ///< Indices of front hit keys - HitSetArray_t<qa::H1D*> fvphHitBackKeyIndex = {nullptr, nullptr}; ///< Indices of back hit keys - - std::vector<HitSetArray_t<qa::H1D*>> fvphHitTime; ///< Time distribution of hits - }; -} // namespace cbm::algo::ca diff --git a/algo/ca/qa/CaOutputQa.cxx b/algo/ca/qa/CaOutputQa.cxx deleted file mode 100644 index ee728fed70..0000000000 --- a/algo/ca/qa/CaOutputQa.cxx +++ /dev/null @@ -1,178 +0,0 @@ -/* Copyright (C) 2023-2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt - SPDX-License-Identifier: GPL-3.0-only - Authors: Sergei Zharko [committer] */ - -/// \file CaOutputQa.cxx -/// \date 20.11.2023 -/// \brief A QA module for CA tracking (implementation) -/// \author S.Zharko <s.zharko@gsi.de> - -#include "CaOutputQa.h" - -#include "CaConstants.h" -#include "CaInputData.h" -#include "CaParameters.h" -#include "CaTrack.h" - -#include <fmt/format.h> - -using cbm::algo::ca::OutputQa; -using cbm::algo::ca::constants::math::Pi; - -// --------------------------------------------------------------------------------------------------------------------- -// -void OutputQa::Init() -{ - using cbm::algo::qa::CanvasConfig; - using cbm::algo::qa::Data; - using cbm::algo::qa::H1D; - using cbm::algo::qa::H2D; - using cbm::algo::qa::PadConfig; - using fmt::format; - - if (!fpSender.get()) { - return; - } - - // ----- Histograms initialization - constexpr int nParPads = 3; - std::array<std::string, knTrkParPoints> vsPointName = {"first", "last"}; - for (int i = 0; i < knTrkParPoints; ++i) { - { - auto sName = format("track_{}_theta", vsPointName[i]); - auto sTitl = format("#theta at {} hit; #theta", vsPointName[i]); - fvphTrkTheta[i] = fQaData.MakeObj<H1D>(sName, sTitl, 62, 0., 90.); - } - { - auto sName = format("track_{}_phi", vsPointName[i]); - auto sTitl = format("#phi at {} hit; #phi", vsPointName[i]); - fvphTrkPhi[i] = fQaData.MakeObj<H1D>(sName, sTitl, 62, -180., 180.); - } - { - auto sName = format("track_{}_thata_phi", vsPointName[i]); - auto sTitl = format("#theta vs #phi at {} hit; #phi; #theta", vsPointName[i]); - fvphTrkPhiTheta[i] = fQaData.MakeObj<H2D>(sName, sTitl, 62, -180., 180., 62, 0., 90.); - } - { - auto sName = format("track_{}_chi2_ndf", vsPointName[i]); - auto sTitl = format("#chi^{{2}}/NDF at {} hit; #chi^{{2}}/NDF", vsPointName[i]); - fvphTrkChi2Ndf[i] = fQaData.MakeObj<H1D>(sName, sTitl, 100, 0., 20.); - } - } - { - int nBins = knStaMax; - double xMin = -0.5; - double xMax = double(knStaMax) - 0.5; - { - auto sName = "track_fst_lst_sta"; - auto sTitl = "First vs. last station index;ID^{last}_{station};ID^{first}_{station}"; - fphTrkFstLstSta = fQaData.MakeObj<H2D>(sName, sTitl, nBins, xMin, xMax, nBins, xMin, xMax); - } - fphTrkNofHits = fQaData.MakeObj<H1D>("n_hits", "Number of hits;N_{hit}", nBins, xMin, xMax); - } - - // ---- Init canvases - { - - // Track parameters at first/last track hit - - for (int i = 0; i < knTrkParPoints; ++i) { - auto sCanvName = fmt::format("ca_trk_{}_hit", vsPointName[i]); - auto sCanvTitl = fmt::format("Track parameters at {} hit", vsPointName[i]); - auto canv = CanvasConfig(sCanvName, sCanvTitl); - std::vector<PadConfig> pads; - pads.reserve(nParPads); - for (int iPad = 0; iPad < nParPads; ++iPad) { - auto& pad = pads.emplace_back(); - pad.SetLog(false, true); - pad.SetGrid(true, true); - } - pads[0].RegisterHistogram(fvphTrkTheta[i], "hist"); - pads[1].RegisterHistogram(fvphTrkPhi[i], "hist"); - pads[2].RegisterHistogram(fvphTrkChi2Ndf[i], "hist"); - for (auto& pad : pads) { - canv.AddPadConfig(pad); - } - fQaData.AddCanvasConfig(canv); - } - - { - auto sCanvName = "ca_trk_theta_phi"; - auto sCanvTitl = "Track #theta vs #phi"; - auto canv = CanvasConfig(sCanvName, sCanvTitl); - for (int i = 0; i < knTrkParPoints; ++i) { - PadConfig pad; - pad.SetLog(false, false, true); - pad.SetGrid(true, true); - pad.RegisterHistogram(fvphTrkPhiTheta[i], "colz"); - canv.AddPadConfig(pad); - } - fQaData.AddCanvasConfig(canv); - } - - // Number of hits in track - { - auto canv = CanvasConfig("ca_nhits_in_trk", "Number of hit in track"); - { - auto pad = PadConfig(); - pad.SetGrid(true, true); - pad.SetLog(false, true); - pad.RegisterHistogram(fphTrkNofHits, "hist"); - canv.AddPadConfig(pad); - } - { - auto pad = PadConfig(); - pad.SetGrid(true, true); - pad.SetLog(false, false, true); - pad.RegisterHistogram(fphTrkFstLstSta, "colz"); - canv.AddPadConfig(pad); - } - fQaData.AddCanvasConfig(canv); - } - } - - fQaData.Init(fpSender); -} - -// --------------------------------------------------------------------------------------------------------------------- -// -void OutputQa::Exec() -{ - if (!fpSender.get()) { - return; - } - - if (!CheckInit()) { - L_(fatal) << "ca::OutputQa: instance is not initialized"; - assert(false); - } - - // Fill track histograms - { - int trkFirstHit = 0; // Index of hit in fpvRecoHits - for (auto trkIt = fpvTracks->begin(); trkIt != fpvTracks->end(); ++trkIt) { - auto& track = *trkIt; - int nHits = track.fNofHits; - // Indices of hits in fpInputData->GetHits() - int iFstHit = (*fpvRecoHits)[trkFirstHit]; - int iLstHit = (*fpvRecoHits)[trkFirstHit + nHits - 1]; - - // Distributions in different track points - for (int ip = 0; ip < knTrkParPoints; ++ip) { - //int iHit = (ip == 0 ? iFstHit : iLstHit); - //const auto& hit = fpInputData->GetHit(iHit); - const auto& trkPar = (ip == 0 ? track.fParFirst : track.fParLast); - fvphTrkTheta[ip]->Fill(trkPar.GetTheta() * 180. / Pi); - fvphTrkPhi[ip]->Fill(trkPar.GetPhi() * 180. / Pi); - fvphTrkPhiTheta[ip]->Fill(trkPar.GetPhi() * 180. / Pi, trkPar.GetTheta() * 180. / Pi); - fvphTrkChi2Ndf[ip]->Fill(trkPar.GetChiSq() / trkPar.GetNdf()); - } - // Other distributions - fphTrkFstLstSta->Fill(fpInputData->GetHit(iLstHit).Station(), fpInputData->GetHit(iFstHit).Station()); - fphTrkNofHits->Fill(nHits); - trkFirstHit += nHits; - } - } - - fQaData.Send(fpSender); -} diff --git a/algo/ca/qa/CaOutputQa.h b/algo/ca/qa/CaOutputQa.h deleted file mode 100644 index 4b089fcf01..0000000000 --- a/algo/ca/qa/CaOutputQa.h +++ /dev/null @@ -1,63 +0,0 @@ -/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt - SPDX-License-Identifier: GPL-3.0-only - Authors: Sergei Zharko [committer] */ - -/// \file CaOutputQa.h -/// \date 20.11.2023 -/// \brief A QA module for CA tracking (header) -/// \author S.Zharko <s.zharko@gsi.de> - -#pragma once - -#include "CaQaBase.h" - -namespace cbm::algo::ca -{ - /// \class cbm::algo::ca::qa::OutputQa - /// \brief OutputQa class for the CA tracking QA (header) - /// - class OutputQa : public QaBase { - public: - /// \brief Default destructor - /// \param pSender Pointer to the histogram sender - OutputQa(std::shared_ptr<HistogramSender> pSender) : QaBase(pSender, "CA/Output"){}; - - /// \brief Constructor from the configuration object - /// \param config QA configuration object - OutputQa() = default; - - /// \brief Copy constructor - OutputQa(const OutputQa&) = delete; - - /// \brief Move constructor - OutputQa(OutputQa&&) = delete; - - /// \brief Destructor - ~OutputQa() = default; - - /// \brief Copy assignment operator - OutputQa& operator=(const OutputQa&) = delete; - - /// \brief Move assignment operator - OutputQa& operator=(OutputQa&&) = delete; - - /// \brief QA execution function - void Exec(); - - /// \brief Initializes the QA - void Init(); - - private: - static constexpr int knTrkParPoints = 2; ///< Number of track points to build par distributions - static constexpr int knStaMax = 16; ///< Max number of stations (histogram binning) - - // Track distributions - std::array<qa::H1D*, knTrkParPoints> fvphTrkTheta = {{0}}; ///< hist: theta at first/last hit - std::array<qa::H1D*, knTrkParPoints> fvphTrkPhi = {{0}}; ///< hist: phi at first/last hit - std::array<qa::H1D*, knTrkParPoints> fvphTrkChi2Ndf = {{0}}; ///< hist: chi2/NDF at first/last hit - std::array<qa::H2D*, knTrkParPoints> fvphTrkPhiTheta = {{0}}; ///< hist: theta vs. phi at first/last hit - - qa::H2D* fphTrkFstLstSta = nullptr; ///< hist: fst vs lst station index - qa::H1D* fphTrkNofHits = nullptr; ///< hist: number of hits in track - }; -} // namespace cbm::algo::ca diff --git a/algo/ca/qa/CaQa.cxx b/algo/ca/qa/CaQa.cxx new file mode 100644 index 0000000000..1dda9fbb17 --- /dev/null +++ b/algo/ca/qa/CaQa.cxx @@ -0,0 +1,402 @@ +/* Copyright (C) 2023-2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt + SPDX-License-Identifier: GPL-3.0-only + Authors: Sergei Zharko [committer] */ + +/// \file CaQa.cxx +/// \date 20.11.2023 +/// \brief A QA module for CA tracking (implementation) +/// \author S.Zharko <s.zharko@gsi.de> + +#include "CaQa.h" + +#include "CaConstants.h" +#include "CaInputData.h" +#include "CaParameters.h" +#include "CaTrack.h" +#include "TrackingDefs.h" + +#include <algorithm> +#include <fstream> +#include <limits> + +#include <fmt/format.h> + +using cbm::algo::ca::Hit; +using cbm::algo::ca::Qa; +using cbm::algo::ca::constants::math::Pi; + +// --------------------------------------------------------------------------------------------------------------------- +// +bool Qa::CheckInit() const +{ + bool res = true; + if (!fpParameters) { + L_(error) << "cbm::algo::ca::OutputQa::Build(): parameters object is undefined"; + res = false; + } + if (!fpInputData) { + L_(error) << "cbm::algo::ca::OutputQa::Build(): input data object is undefined"; + res = false; + } + if (!fpvTracks) { + L_(error) << "cbm::algo::ca::OutputQa::Build(): track vector is undefined"; + res = false; + } + if (!fpvRecoHits) { + L_(error) << "cbm::algo::ca::OutputQa::Build(): used hit index vector is undefined"; + res = false; + } + return res; +} + +// --------------------------------------------------------------------------------------------------------------------- +// +void Qa::Init() +{ + using cbm::algo::qa::CanvasConfig; + using cbm::algo::qa::Data; + using cbm::algo::qa::H1D; + using cbm::algo::qa::H2D; + using cbm::algo::qa::PadConfig; + using cbm::algo::qa::Prof2D; + using fmt::format; + + if (!fpSender.get()) { + return; + } + + // ******************** + // ** Hit distributions + int nSt = fpParameters->GetNstationsActive(); // Number of active tracking stations + { + // Occupancy distributions + fvphHitOccupXY.resize(nSt); + fvphHitOccupZX.resize(nSt); + if (kDebug) { + fvphHitOccupZY.resize(nSt); + fvphHitUsageXY.resize(nSt); + } + + // Station sizes in transverse plane + std::vector<double> vMinX(nSt); + std::vector<double> vMaxX(nSt); + std::vector<double> vMinY(nSt); + std::vector<double> vMaxY(nSt); + + int nBinsXY = 200; + int nBinsZ = 400; + + for (int iSt = 0; iSt < nSt; ++iSt) { + const auto& station = fpParameters->GetStation(iSt); + vMaxX[iSt] = station.GetXmax<double>(); + vMinX[iSt] = station.GetXmin<double>(); + vMaxY[iSt] = station.GetYmax<double>(); + vMinY[iSt] = station.GetYmin<double>(); + double dy = (vMaxY[iSt] - vMinY[iSt]) * kXYZMargin; + double dx = (vMaxX[iSt] - vMinX[iSt]) * kXYZMargin; + vMaxX[iSt] += dx; + vMinX[iSt] -= dx; + vMaxY[iSt] += dy; + vMinY[iSt] -= dy; + } + // Station max + double xMinA = *std::min_element(vMinX.begin(), vMinX.end()); + double xMaxA = *std::max_element(vMaxX.begin(), vMaxX.end()); + double yMinA = *std::min_element(vMinY.begin(), vMinY.end()); + double yMaxA = *std::max_element(vMaxY.begin(), vMaxY.end()); + double zMinA = fpParameters->GetStation(0).GetZ<double>(); + double zMaxA = fpParameters->GetStation(nSt - 1).GetZ<double>(); + { + double dz = (zMaxA - zMinA) * kXYZMargin; + zMinA -= dz; + zMaxA += dz; + } + for (int iSt = 0; iSt < nSt; ++iSt) { + int iStGeo = fpParameters->GetActiveToGeoMap()[iSt]; + auto [detID, iStLoc] = fpParameters->GetGeoToLocalIdMap()[iStGeo]; + for (auto hitSet : kHitSets) { + auto setNm = EHitSet::Input == hitSet ? "input" : "used"; + auto setTl = EHitSet::Input == hitSet ? "Input" : "Used"; + { + auto name = format("hit_{}_occup_xy_sta_{}", setNm, iSt); + auto titl = format("{} hit occupancy in XY plane for station {} ({}{});x [cm];y [cm]", setTl, iSt, + kDetName[detID], iStLoc); + fvphHitOccupXY[iSt][hitSet] = + fQaData.MakeObj<H2D>(name, titl, nBinsXY, vMinX[iSt], vMaxX[iSt], nBinsXY, vMinY[iSt], vMaxY[iSt]); + } + { + auto name = format("hit_{}_occup_zx_sta_{}", setNm, iSt); + auto titl = format("{} hit occupancy in ZX plane for station {} ({}{});z [cm];x [cm]", setTl, iSt, + kDetName[detID], iStLoc); + fvphHitOccupZX[iSt][hitSet] = fQaData.MakeObj<H2D>(name, titl, nBinsZ, zMinA, zMaxA, nBinsXY, xMinA, xMaxA); + } + if (kDebug) { + auto name = format("hit_{}_occup_zy_sta_{}", setNm, iSt); + auto titl = format("{} hit occupancy in ZY plane for station {} ({}{});z [cm];y [cm]", setTl, iSt, + kDetName[detID], iStLoc); + fvphHitOccupZY[iSt][hitSet] = fQaData.MakeObj<H2D>(name, titl, nBinsZ, zMinA, zMaxA, nBinsXY, yMinA, yMaxA); + } + } + if (kDebug) { + auto name = format("hit_usage_xy_sta_{}", iSt); + auto titl = format("Hit usage in XY plane for station {} ({}{});x [cm];y [cm]", iSt, kDetName[detID], iStLoc); + fvphHitUsageXY[iSt] = + fQaData.MakeObj<Prof2D>(name, titl, nBinsXY, vMinX[iSt], vMaxX[iSt], nBinsXY, vMinY[iSt], vMaxY[iSt], 0., 1.); + } + } + if (kDebug) { + for (auto hitSet : kHitSets) { + constexpr int NBins = 1000000; + auto setNm = EHitSet::Input == hitSet ? "input" : "used"; + auto setTl = EHitSet::Input == hitSet ? "Input" : "Used"; + { + auto name = format("hit_{}_front_key_index", setNm); + auto titl = format("{} hit front key index;ID_{{key}}/N_{{keys}};Count", setTl); + fvphHitFrontKeyIndex[hitSet] = fQaData.MakeObj<H1D>(name, titl, NBins, 0., 1.); + } + { + auto name = format("hit_{}_back_key_index", setNm); + auto titl = format("{} hit back key index;ID_{{key}}/N_{{keys}};Count", setTl); + fvphHitBackKeyIndex[hitSet] = fQaData.MakeObj<H1D>(name, titl, NBins, 0., 1.); + } + } + } + + // Hit time distributions + if (kDebug) { + fvphHitTime.resize(nSt + 1); // [nSt] - total over stations + for (int iSt = 0; iSt < nSt + 1; ++iSt) { + auto staNm = iSt == nSt ? "" : format("_sta_{}", iSt); + auto staTl = iSt == nSt ? "" : format(" in station {}", iSt); + for (auto hitSet : kHitSets) { + auto setNm = EHitSet::Input == hitSet ? "input" : "used"; + auto setTl = EHitSet::Input == hitSet ? "Input" : "Used"; + auto name = format("hit_{}_rel_time{}", setNm, staNm); + auto titl = format("{} hit relative time{}; #delta t_{{hit}};Count", setTl, staTl); + fvphHitTime[iSt][hitSet] = fQaData.MakeObj<H1D>(name, titl, 10000, -0.1, 1.1); + } + } + } + } + + // ********************** + // ** Track distributions + std::vector<std::string> vsPointName = {"first", "last"}; + for (int i = 0; i < knTrkParPoints; ++i) { + if (!kDebug && i > 0) { + continue; + } + { + auto sName = format("track_{}_theta", vsPointName[i]); + auto sTitl = format("#theta at {} hit; #theta", vsPointName[i]); + fvphTrkTheta[i] = fQaData.MakeObj<H1D>(sName, sTitl, 62, 0., 90.); + } + { + auto sName = format("track_{}_phi", vsPointName[i]); + auto sTitl = format("#phi at {} hit; #phi", vsPointName[i]); + fvphTrkPhi[i] = fQaData.MakeObj<H1D>(sName, sTitl, 62, -180., 180.); + } + { + auto sName = format("track_{}_thata_phi", vsPointName[i]); + auto sTitl = format("#theta vs #phi at {} hit; #phi; #theta", vsPointName[i]); + fvphTrkPhiTheta[i] = fQaData.MakeObj<H2D>(sName, sTitl, 62, -180., 180., 62, 0., 90.); + } + { + auto sName = format("track_{}_chi2_ndf", vsPointName[i]); + auto sTitl = format("#chi^{{2}}/NDF at {} hit; #chi^{{2}}/NDF", vsPointName[i]); + fvphTrkChi2Ndf[i] = fQaData.MakeObj<H1D>(sName, sTitl, 100, 0., 20.); + } + } + { + int nBins = knStaMax; + double xMin = -0.5; + double xMax = double(knStaMax) - 0.5; + { + auto sName = "track_fst_lst_sta"; + auto sTitl = "First vs. last station index;ID^{last}_{station};ID^{first}_{station}"; + fphTrkFstLstSta = fQaData.MakeObj<H2D>(sName, sTitl, nBins, xMin, xMax, nBins, xMin, xMax); + } + fphTrkNofHits = fQaData.MakeObj<H1D>("n_hits", "Number of hits;N_{hit}", nBins, xMin, xMax); + } + + // ---- Init canvases + { + // Input hits canvas + { + for (auto hitSet : kHitSets) { + auto setNm = EHitSet::Input == hitSet ? "input" : "used"; + auto setTl = EHitSet::Input == hitSet ? "Input" : "Used"; + { // XY + auto name = format("ca_hit_{}_occupancy_xy", setNm); + auto titl = format("{} hit occupancy in different stations in XY plane", setNm); + auto canv = CanvasConfig(name, titl); + for (int iSt = 0; iSt < nSt; ++iSt) { + auto pad = PadConfig(); + pad.RegisterHistogram(fvphHitOccupXY[iSt][hitSet], "colz"); + canv.AddPadConfig(pad); + } + fQaData.AddCanvasConfig(canv); + } + { // ZX and ZY + auto name = format("ca_hit_{}_occupancy_zx_zy", setNm); + auto titl = format("{} hit occupancy in different stations in ZX and ZY planes", setTl); + auto canv = CanvasConfig(name, titl); + { // ZX + auto pad = PadConfig(); + for (int iSt = 0; iSt < nSt; ++iSt) { + pad.RegisterHistogram(fvphHitOccupZX[iSt][hitSet], (iSt == 0 ? "colz" : "cols same")); + } + canv.AddPadConfig(pad); + } + if (kDebug) { // ZY + auto pad = PadConfig(); + for (int iSt = 0; iSt < nSt; ++iSt) { + pad.RegisterHistogram(fvphHitOccupZY[iSt][hitSet], (iSt == 0 ? "colz" : "cols same")); + } + canv.AddPadConfig(pad); + } + fQaData.AddCanvasConfig(canv); + } + } + if (kDebug) { + auto name = format("ca_hit_usage_xy"); + auto titl = format("Hit usage in different stations in XY plane"); + auto canv = CanvasConfig(name, titl); + for (int iSt = 0; iSt < nSt; ++iSt) { + auto pad = PadConfig(); + pad.RegisterHistogram(fvphHitUsageXY[iSt], "colz"); + canv.AddPadConfig(pad); + } + fQaData.AddCanvasConfig(canv); + } + } + + // Tracks canvas + { + auto canv = CanvasConfig("ca_tracks", "Tracking output QA"); + { + auto pad = PadConfig(true, true, false, false, true); + pad.RegisterHistogram(fvphTrkPhiTheta[0], "colz"); + canv.AddPadConfig(pad); + } + { + auto pad = PadConfig(true, true, false, true, false); + pad.RegisterHistogram(fvphTrkChi2Ndf[0], "hist"); + canv.AddPadConfig(pad); + } + { + auto pad = PadConfig(true, true, false, true, false); + pad.RegisterHistogram(fphTrkNofHits, "hist"); + canv.AddPadConfig(pad); + } + { + auto pad = PadConfig(true, true, false, false, true); + pad.RegisterHistogram(fphTrkFstLstSta, "colz"); + canv.AddPadConfig(pad); + } + fQaData.AddCanvasConfig(canv); + } + fQaData.Init(fpSender); + } +} + +// --------------------------------------------------------------------------------------------------------------------- +// +void Qa::Exec() +{ + if (!fpSender.get()) { + return; + } + + if (!CheckInit()) { + L_(fatal) << "ca::Qa: instance is not initialized"; + assert(false); + } + + int nHitsInput = fpInputData->GetNhits(); + // Map: if hit used in tracking + std::vector<unsigned char> vbHitUsed(nHitsInput, false); + for (int iH : (*fpvRecoHits)) { + vbHitUsed[iH] = true; + } + + // Calculate max and min hit time + if constexpr (kDebug) { + const auto& hits = fpInputData->GetHits(); + auto [minTimeIt, maxTimeIt] = + std::minmax_element(hits.cbegin(), hits.cend(), [](const auto& h1, const auto& h2) { return h1.T() < h2.T(); }); + fMinHitTime = minTimeIt->T(); + fMaxHitTime = maxTimeIt->T(); + } + + // Fill input hit histograms + { + for (int iH = 0; iH < nHitsInput; ++iH) { + const auto& hit = fpInputData->GetHit(iH); + FillHitDistributionsForHitSet(EHitSet::Input, hit); + if (vbHitUsed[iH]) { + FillHitDistributionsForHitSet(EHitSet::Used, hit); + } + if constexpr (kDebug) { + int iSt = hit.Station(); + double x = hit.X(); + double y = hit.Y(); + fvphHitUsageXY[iSt]->Fill(x, y, static_cast<double>(vbHitUsed[iH])); + } + } + } + + // Fill track histograms + { + int trkFirstHit = 0; // Index of hit in fpvRecoHits + for (auto trkIt = fpvTracks->begin(); trkIt != fpvTracks->end(); ++trkIt) { + auto& track = *trkIt; + int nHits = track.fNofHits; + // Indices of hits in fpInputData->GetHits() + int iFstHit = (*fpvRecoHits)[trkFirstHit]; + int iLstHit = (*fpvRecoHits)[trkFirstHit + nHits - 1]; + + // Distributions in different track points + for (int ip = 0; ip < knTrkParPoints; ++ip) { + if (!kDebug && ip > 0) { + continue; + } + //int iHit = (ip == 0 ? iFstHit : iLstHit); + //const auto& hit = fpInputData->GetHit(iHit); + const auto& trkPar = (ip == 0 ? track.fParFirst : track.fParLast); + fvphTrkTheta[ip]->Fill(trkPar.GetTheta() * 180. / Pi); + fvphTrkPhi[ip]->Fill(trkPar.GetPhi() * 180. / Pi); + fvphTrkPhiTheta[ip]->Fill(trkPar.GetPhi() * 180. / Pi, trkPar.GetTheta() * 180. / Pi); + fvphTrkChi2Ndf[ip]->Fill(trkPar.GetChiSq() / trkPar.GetNdf()); + } + // Other distributions + fphTrkFstLstSta->Fill(fpInputData->GetHit(iLstHit).Station(), fpInputData->GetHit(iFstHit).Station()); + fphTrkNofHits->Fill(nHits); + trkFirstHit += nHits; + } + } + + fQaData.Send(fpSender); +} + +// --------------------------------------------------------------------------------------------------------------------- +// +void Qa::FillHitDistributionsForHitSet(Qa::EHitSet hitSet, const Hit& hit) +{ + int nSt = fpParameters->GetNstationsActive(); + int iSt = hit.Station(); + double x = hit.X(); + double y = hit.Y(); + double z = hit.Z(); + fvphHitOccupXY[iSt][hitSet]->Fill(x, y); + fvphHitOccupZX[iSt][hitSet]->Fill(z, x); + if constexpr (kDebug) { + fvphHitOccupZY[iSt][hitSet]->Fill(z, y); + auto nKeys = static_cast<double>(fpInputData->GetNhitKeys()); + fvphHitFrontKeyIndex[hitSet]->Fill(hit.FrontKey() / nKeys); + fvphHitBackKeyIndex[hitSet]->Fill(hit.BackKey() / nKeys); + double relTime = (hit.T() - fMinHitTime) / (fMaxHitTime - fMinHitTime); + fvphHitTime[iSt][hitSet]->Fill(relTime); + fvphHitTime[nSt][hitSet]->Fill(relTime); + } +} diff --git a/algo/ca/qa/CaQa.h b/algo/ca/qa/CaQa.h new file mode 100644 index 0000000000..dc687a4221 --- /dev/null +++ b/algo/ca/qa/CaQa.h @@ -0,0 +1,156 @@ +/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt + SPDX-License-Identifier: GPL-3.0-only + Authors: Sergei Zharko [committer] */ + +/// \file CaQa.h +/// \date 20.11.2023 +/// \brief A QA module for CA tracking (header) +/// \author S.Zharko <s.zharko@gsi.de> + +#pragma once + +#include "CaEnumArray.h" +#include "CaHit.h" // for HitIndex_t +#include "CaTimesliceHeader.h" +#include "CaVector.h" +#include "QaData.h" // QA data + +namespace cbm::algo +{ + namespace qa + { + class H1D; + class H2D; + } // namespace qa + namespace ca + { + template<typename DataT> + class Parameters; + class InputData; + class Track; + } // namespace ca +} // namespace cbm::algo + +namespace cbm::algo::ca +{ + /// \class cbm::algo::ca::qa::Qa + /// \brief Qa class for the CA tracking QA (header) + /// + class Qa { + /// \brief Hit set entries + enum class EHitSet + { + Input, ///< Input hits + Used, ///< Hits used in tracks + kEND + }; + + /// \brief Definition of enum array over EHitSet entries + template<typename T> + using HitSetArray_t = EnumArray<EHitSet, T>; + + /// \brief Array of EHitSet entries for iteration + static constexpr HitSetArray_t<EHitSet> kHitSets = {EHitSet::Input, EHitSet::Used}; + + public: + /// \brief Default destructor + /// \param pSender Pointer to the histogram sender + Qa(std::shared_ptr<HistogramSender> pSender) : fQaData("CA"), fpSender(pSender) {} + + /// \brief Constructor from the configuration object + /// \param config QA configuration object + Qa() = default; + + /// \brief Copy constructor + Qa(const Qa&) = delete; + + /// \brief Move constructor + Qa(Qa&&) = delete; + + /// \brief Destructor + ~Qa() = default; + + /// \brief Copy assignment operator + Qa& operator=(const Qa&) = delete; + + /// \brief Move assignment operator + Qa& operator=(Qa&&) = delete; + + /// \brief QA execution function + void Exec(); + + /// \brief Initializes the QA + void Init(); + + /// \brief Check initialization + /// \return true All variables are initialized + /// \return false Some of are not initialized + bool CheckInit() const; + + /// \brief Checks, if the histogram sender is defined + bool IsSenderDefined() const { return static_cast<bool>(fpSender.get()); } + + /// \brief Registers tracking input data object + /// \note Call per TS + void RegisterInputData(const InputData* pInputData) { fpInputData = pInputData; } + + /// \brief Registers track vector + /// \note Call per TS + void RegisterTracks(const Vector<Track>* pvTracks) { fpvTracks = pvTracks; } + + /// \brief Registers reco hits indices vector + /// \note Call per TS + void RegisterRecoHitIndices(const Vector<HitIndex_t>* pvRecoHits) { fpvRecoHits = pvRecoHits; } + + /// \brief Registers tracking parameters object + /// \note Call per run + void RegisterParameters(const Parameters<fvec>* pParameters) { fpParameters = pParameters; } + + private: + /// \brief Fills hit distributions + /// \param hitSet Hit set enum entry + /// \param hit Reference to hit + void FillHitDistributionsForHitSet(EHitSet hitSet, const ca::Hit& hit); + + // parameters + static constexpr double kXYZMargin = 0.05; ///< Margin for occupancy distributions in XY plane + static constexpr int knHitSets = 2; ///< Number of hit sets: input/used + static constexpr int knTrkParPoints = 2; ///< Number of track points to build par distributions + static constexpr int knStaMax = 16; ///< Max number of stations (histogram binning) + static constexpr bool kDebug = false; ///< Additional histograms + + double fMinHitTime = std::numeric_limits<double>::max(); + double fMaxHitTime = std::numeric_limits<double>::lowest(); + + // utility + qa::Data fQaData; ///< QA data + + std::shared_ptr<HistogramSender> fpSender = nullptr; ///< Histogram sender + const Parameters<fvec>* fpParameters = nullptr; ///< Pointer to tracking parameters + const InputData* fpInputData = nullptr; ///< Pointer to input data + const Vector<Track>* fpvTracks = nullptr; ///< Pointer to tracks vector + const Vector<HitIndex_t>* fpvRecoHits = nullptr; ///< Pointer to reco hit indices + + // Hit distributions + using OccupHistContainer_t = std::vector<HitSetArray_t<qa::H2D*>>; + OccupHistContainer_t fvphHitOccupXY; ///< hist: Hit occupancy in different stations in XY plane + OccupHistContainer_t fvphHitOccupZX; ///< hist: Hit occupancy in different stations in ZX plane + OccupHistContainer_t fvphHitOccupZY; ///< hist: Hit occupancy in different stations in ZY plane + + std::vector<qa::Prof2D*> fvphHitUsageXY; ///< prof: Hit usage in different stations in XY plane + + HitSetArray_t<qa::H1D*> fvphHitFrontKeyIndex = {nullptr, nullptr}; ///< Indices of front hit keys + HitSetArray_t<qa::H1D*> fvphHitBackKeyIndex = {nullptr, nullptr}; ///< Indices of back hit keys + + std::vector<HitSetArray_t<qa::H1D*>> fvphHitTime; ///< Time distribution of hits + + // Track distributions + std::array<qa::H1D*, knTrkParPoints> fvphTrkTheta = {{0}}; ///< hist: theta at first/last hit + std::array<qa::H1D*, knTrkParPoints> fvphTrkPhi = {{0}}; ///< hist: phi at first/last hit + std::array<qa::H1D*, knTrkParPoints> fvphTrkChi2Ndf = {{0}}; ///< hist: chi2/NDF at first/last hit + std::array<qa::H2D*, knTrkParPoints> fvphTrkPhiTheta = {{0}}; ///< hist: theta vs. phi at first/last hit + + qa::H2D* fphTrkFstLstSta = nullptr; ///< hist: fst vs lst station index + qa::H1D* fphTrkNofHits = nullptr; ///< hist: number of hits in track + }; +} // namespace cbm::algo::ca diff --git a/algo/ca/qa/CaQaBase.cxx b/algo/ca/qa/CaQaBase.cxx deleted file mode 100644 index da2810670e..0000000000 --- a/algo/ca/qa/CaQaBase.cxx +++ /dev/null @@ -1,36 +0,0 @@ -/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt - SPDX-License-Identifier: GPL-3.0-only - Authors: Sergei Zharko [committer] */ - -/// \file CaQaBase.h -/// \date 01.03.2024 -/// \brief Base class for tracking QA (source) -/// \author Sergei Zharko <s.zharko@gsi.de> - -#include "CaQaBase.h" - -using cbm::algo::ca::QaBase; - -// --------------------------------------------------------------------------------------------------------------------- -// -bool QaBase::CheckInit() const -{ - bool res = true; - if (!fpParameters) { - L_(error) << "cbm::algo::ca::OutputQa::Build(): parameters object is undefined"; - res = false; - } - if (!fpInputData) { - L_(error) << "cbm::algo::ca::OutputQa::Build(): input data object is undefined"; - res = false; - } - if (!fpvTracks) { - L_(error) << "cbm::algo::ca::OutputQa::Build(): track vector is undefined"; - res = false; - } - if (!fpvRecoHits) { - L_(error) << "cbm::algo::ca::OutputQa::Build(): used hit index vector is undefined"; - res = false; - } - return res; -} diff --git a/algo/ca/qa/CaQaBase.h b/algo/ca/qa/CaQaBase.h deleted file mode 100644 index e1fbeb6d1c..0000000000 --- a/algo/ca/qa/CaQaBase.h +++ /dev/null @@ -1,100 +0,0 @@ -/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt - SPDX-License-Identifier: GPL-3.0-only - Authors: Sergei Zharko [committer] */ - -/// \file CaQaBase.h -/// \date 29.02.2024 -/// \brief Base class for tracking QA (header) -/// \author Sergei Zharko <s.zharko@gsi.de> - -#pragma once - -#include "CaHit.h" // for HitIndex_t -#include "CaTimesliceHeader.h" -#include "CaVector.h" -#include "QaData.h" // QA data - -namespace cbm::algo -{ - namespace qa - { - class H1D; - class H2D; - } // namespace qa - namespace ca - { - template<typename DataT> - class Parameters; - class InputData; - class Track; - } // namespace ca -} // namespace cbm::algo - -namespace cbm::algo::ca -{ - /// \class cbm::algo::ca::qa::QaBase - /// \brief Base class for CA tracking QA - /// - /// The class contains set of pointers to data structures, which are needed to build the histograms. - class QaBase { - public: - /// \brief Default destructor - /// \param pSender Pointer to the histogram sender - /// \param dirname Directory name in the output file - QaBase(std::shared_ptr<HistogramSender> pSender, const std::string& dirname) : fQaData(dirname), fpSender(pSender) - { - } - - /// \brief Constructor from the configuration object - /// \param config QA configuration object - QaBase() = default; - - /// \brief Copy constructor - QaBase(const QaBase&) = delete; - - /// \brief Move constructor - QaBase(QaBase&&) = delete; - - /// \brief Destructor - ~QaBase() = default; - - /// \brief Copy assignment operator - QaBase& operator=(const QaBase&) = delete; - - /// \brief Move assignment operator - QaBase& operator=(QaBase&&) = delete; - - /// \brief Check initialization - /// \return true All variables are initialized - /// \return false Some of are not initialized - bool CheckInit() const; - - /// \brief Checks, if the histogram sender is defined - bool IsSenderDefined() const { return static_cast<bool>(fpSender.get()); } - - /// \brief Registers tracking input data object - /// \note Call per TS - void RegisterInputData(const InputData* pInputData) { fpInputData = pInputData; } - - /// \brief Registers track vector - /// \note Call per TS - void RegisterTracks(const Vector<Track>* pvTracks) { fpvTracks = pvTracks; } - - /// \brief Registers reco hits indices vector - /// \note Call per TS - void RegisterRecoHitIndices(const Vector<HitIndex_t>* pvRecoHits) { fpvRecoHits = pvRecoHits; } - - /// \brief Registers tracking parameters object - /// \note Call per run - void RegisterParameters(const Parameters<fvec>* pParameters) { fpParameters = pParameters; } - - protected: - qa::Data fQaData; ///< QA data - - std::shared_ptr<HistogramSender> fpSender = nullptr; ///< Histogram sender - const Parameters<fvec>* fpParameters = nullptr; ///< Pointer to tracking parameters - const InputData* fpInputData = nullptr; ///< Pointer to input data - const Vector<Track>* fpvTracks = nullptr; ///< Pointer to tracks vector - const Vector<HitIndex_t>* fpvRecoHits = nullptr; ///< Pointer to reco hit indices - }; -} // namespace cbm::algo::ca diff --git a/algo/qa/PadConfig.h b/algo/qa/PadConfig.h index 1c4d3cb992..a7dc7469d4 100644 --- a/algo/qa/PadConfig.h +++ b/algo/qa/PadConfig.h @@ -28,6 +28,22 @@ namespace cbm::algo::qa /// \brief Constructor PadConfig() = default; + /// \brief Constructor from parameters + /// \param gridX + /// \param gridY + /// \param logX + /// \param logY + /// \param logZ + PadConfig(bool gridX, bool gridY, bool logX, bool logY, bool logZ) + : fbGridX(gridX) + , fbGridY(gridY) + , fbLogX(logX) + , fbLogY(logY) + , fbLogZ(logZ) + { + } + + /// \brief Copy constructor PadConfig(const PadConfig&) = default; diff --git a/core/base/utils/CbmMcbmUtils.cxx b/core/base/utils/CbmMcbmUtils.cxx index 732256777e..9fc673f146 100644 --- a/core/base/utils/CbmMcbmUtils.cxx +++ b/core/base/utils/CbmMcbmUtils.cxx @@ -59,8 +59,7 @@ namespace cbm } else if (2918 <= ulRunId) { /// Dummy needed to run the unpack macro until we have a setup ready - /// FIXME: replace with 2024 nickel setup!!!!! - sSetupName = "mcbm_beam_2024_03_22_gold"; + sSetupName = "mcbm_beam_2024_05_08_nickel"; } else { /// Missing runs, exception there to force implementation and support from users side -- GitLab