From b409eafebbf73273a942d45558d3b910d40e88c5 Mon Sep 17 00:00:00 2001 From: "s.zharko@gsi.de" <s.zharko@gsi.de> Date: Mon, 28 Aug 2023 13:20:15 +0200 Subject: [PATCH] CA: QA task for setup in tracking --- macro/L1/configs/ca_params_global.yaml | 4 +- macro/L1/configs/ca_params_mcbm.yaml | 4 +- macro/L1/configs/ca_params_sts.yaml | 5 +- macro/mcbm/mcbm_qa.C | 8 ++ reco/L1/CMakeLists.txt | 1 + reco/L1/CbmL1DetectorID.h | 2 + reco/L1/qa/CbmCaInputQaSetup.cxx | 174 +++++++++++++++++++++++++ reco/L1/qa/CbmCaInputQaSetup.h | 145 +++++++++++++++++++++ 8 files changed, 336 insertions(+), 7 deletions(-) create mode 100644 reco/L1/qa/CbmCaInputQaSetup.cxx create mode 100644 reco/L1/qa/CbmCaInputQaSetup.h diff --git a/macro/L1/configs/ca_params_global.yaml b/macro/L1/configs/ca_params_global.yaml index d633f7ce83..faf1ca0126 100644 --- a/macro/L1/configs/ca_params_global.yaml +++ b/macro/L1/configs/ca_params_global.yaml @@ -20,9 +20,9 @@ ca: # interface classes). If there is no index provided, the whole detector subsystem is skept. # Examples: # 1) Turn the first and the second STS in the geometry and the full TRD detector off - # inactive_stations: [STS0', 'STS1', 'TRD'] + # inactive_stations: [STS:0', 'STS:1', 'TRD'] # 2) Turn first TOF station in the geometry off - # inactive_stations: [TOF0] + # inactive_stations: [TOF:0] # 3) Turn the full TOF off # inactive_stations: [TOF] inactive_stations: [] diff --git a/macro/L1/configs/ca_params_mcbm.yaml b/macro/L1/configs/ca_params_mcbm.yaml index 4593b81124..260f7219e8 100644 --- a/macro/L1/configs/ca_params_mcbm.yaml +++ b/macro/L1/configs/ca_params_mcbm.yaml @@ -20,9 +20,9 @@ ca: # interface classes). If there is no index provided, the whole detector subsystem is skept. # Examples: # 1) Turn the first and the second STS in the geometry and the full TRD detector off - # inactive_stations: [STS0', 'STS1', 'TRD'] + # inactive_stations: ['STS:0', 'STS:1', 'TRD'] # 2) Turn first TOF station in the geometry off - # inactive_stations: [TOF0] + # inactive_stations: [TOF:0] # 3) Turn the full TOF off # inactive_stations: [TOF] #inactive_stations: ['MUCH'] diff --git a/macro/L1/configs/ca_params_sts.yaml b/macro/L1/configs/ca_params_sts.yaml index 4172e54cbf..7350f5cd7e 100644 --- a/macro/L1/configs/ca_params_sts.yaml +++ b/macro/L1/configs/ca_params_sts.yaml @@ -20,12 +20,11 @@ ca: # interface classes). If there is no index provided, the whole detector subsystem is skept. # Examples: # 1) Turn the first and the second STS in the geometry and the full TRD detector off - # inactive_stations: [STS0', 'STS1', 'TRD'] + # inactive_stations: [STS:0', 'STS:1', 'TRD'] # 2) Turn first TOF station in the geometry off - # inactive_stations: [TOF0] + # inactive_stations: [TOF:0] # 3) Turn the full TOF off # inactive_stations: [TOF] - #inactive_stations: ['MUCH'] inactive_stations: [] # Random seed diff --git a/macro/mcbm/mcbm_qa.C b/macro/mcbm/mcbm_qa.C index c048add509..d4ec8fcedd 100644 --- a/macro/mcbm/mcbm_qa.C +++ b/macro/mcbm/mcbm_qa.C @@ -246,6 +246,14 @@ void mcbm_qa(Int_t nEvents = 0, TString dataset = "data/mcbm_beam_2020_03_test", TString caParFile = recFile; caParFile.ReplaceAll(".root", ".L1Parameters.dat"); + auto* pCaInputQaSetup = new cbm::ca::InputQaSetup(verbose, bUseMC); + pCaInputQaSetup->SetDetectorFlag(L1DetectorID::kSts, bUseSts); + pCaInputQaSetup->SetDetectorFlag(L1DetectorID::kMuch, bUseMuch); + pCaInputQaSetup->SetDetectorFlag(L1DetectorID::kTrd, bUseTrd); + pCaInputQaSetup->SetDetectorFlag(L1DetectorID::kTof, bUseTof); + pCaInputQaSetup->ReadParameters(caParFile.Data()); + run->AddTask(pCaInputQaSetup); + auto* pCaOutputQa = new cbm::ca::OutputQa(verbose, bUseMC); pCaOutputQa->SetMcbmTrackingMode(); pCaOutputQa->ReadParameters(caParFile.Data()); diff --git a/reco/L1/CMakeLists.txt b/reco/L1/CMakeLists.txt index 080760aca9..c5e6357083 100644 --- a/reco/L1/CMakeLists.txt +++ b/reco/L1/CMakeLists.txt @@ -84,6 +84,7 @@ set(SRCS qa/CbmCaInputQaMuch.cxx qa/CbmCaInputQaTrd.cxx qa/CbmCaInputQaTof.cxx + qa/CbmCaInputQaSetup.cxx qa/CbmCaOutputQa.cxx qa/CbmCaTrackTypeQa.cxx qa/CbmCaTrackFitQa.cxx diff --git a/reco/L1/CbmL1DetectorID.h b/reco/L1/CbmL1DetectorID.h index 7c21822acf..6890d9bd08 100644 --- a/reco/L1/CbmL1DetectorID.h +++ b/reco/L1/CbmL1DetectorID.h @@ -54,6 +54,8 @@ namespace cbm::ca template<typename T> using DetIdArr_t = L1EArray<L1DetectorID, T>; + /// @brief List of + /// @struct DetIdTypeArr_t /// @brief Array of types, indexed by L1DetectorID enum /// diff --git a/reco/L1/qa/CbmCaInputQaSetup.cxx b/reco/L1/qa/CbmCaInputQaSetup.cxx new file mode 100644 index 0000000000..f2f8a1844a --- /dev/null +++ b/reco/L1/qa/CbmCaInputQaSetup.cxx @@ -0,0 +1,174 @@ +/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt + SPDX-License-Identifier: GPL-3.0-only + Authors: Sergei Zharko [committer] */ + +/// @file CbmCaInputQaSetup.cxx +/// @brief QA task for tracking detector interfaces (implementation) +/// @since 28.08.2023 +/// @author S.Zharko <s.zharko@gsi.de> + +#include "CbmCaInputQaSetup.h" + +#include "CbmMCDataManager.h" + +#include "FairRootManager.h" + +#include "L1InitManager.h" + +using cbm::ca::InputQaSetup; + +// --------------------------------------------------------------------------------------------------------------------- +// +InputQaSetup::InputQaSetup(int verbose, bool isMCUsed) : CbmQaTask("CbmCaInputQaSetup", verbose, isMCUsed) {} + +// --------------------------------------------------------------------------------------------------------------------- +// +bool InputQaSetup::Check() { return true; } + + +// --------------------------------------------------------------------------------------------------------------------- +// +void InputQaSetup::CheckInit() const +{ + if (IsMCUsed() && !fpMCEventHeader) { throw std::logic_error("MC event header branch is unavailable"); } + for (int iD = 0; iD < static_cast<int>(L1DetectorID::kEND); ++iD) { + if (fvbUseDet[iD]) { + if (!fvpBrHits[iD]) { throw std::logic_error(Form("Hit branch is unavailable for %s", kDetName[iD])); } + if (IsMCUsed() && !fvpBrPoints[iD]) { + throw std::logic_error(Form("MC point branch is unavailable for %s", kDetName[iD])); + } + } + } +} + +// --------------------------------------------------------------------------------------------------------------------- +// +void InputQaSetup::FillHistograms() +{ + if (fvbUseDet[L1DetectorID::kMvd]) { this->FillHistogramsDet<L1DetectorID::kMvd>(); } + if (fvbUseDet[L1DetectorID::kSts]) { this->FillHistogramsDet<L1DetectorID::kSts>(); } + if (fvbUseDet[L1DetectorID::kMuch]) { this->FillHistogramsDet<L1DetectorID::kMuch>(); } + if (fvbUseDet[L1DetectorID::kTrd]) { this->FillHistogramsDet<L1DetectorID::kTrd>(); } + if (fvbUseDet[L1DetectorID::kTof]) { this->FillHistogramsDet<L1DetectorID::kTof>(); } +} + +// --------------------------------------------------------------------------------------------------------------------- +// +InitStatus InputQaSetup::InitDataBranches() +try { + LOG(info) << fName << ": initializing... "; + + auto pFairManager = FairRootManager::Instance(); + assert(pFairManager); + + auto pMcManager = IsMCUsed() ? dynamic_cast<CbmMCDataManager*>(pFairManager->GetObject("MCDataManager")) : nullptr; + fpMCEventHeader = IsMCUsed() ? pMcManager->GetObject("MCEventHeader.") : nullptr; + + fvpBrPoints.fill(nullptr); + fvpBrHits.fill(nullptr); + + auto InitDetInterface = [&](CbmTrackingDetectorInterfaceBase* ifs, L1DetectorID detID) { + if (fvbUseDet[detID]) { fvpDetInterface[detID] = ifs; } + }; + + auto InitHitBranch = [&](const char* brName, L1DetectorID detID) { + if (fvbUseDet[detID]) { fvpBrHits[detID] = dynamic_cast<TClonesArray*>(pFairManager->GetObject(brName)); } + }; + + auto InitPointBranch = [&](const char* brName, L1DetectorID detID) { + if (IsMCUsed() && fvbUseDet[detID]) { fvpBrPoints[detID] = pMcManager->InitBranch(brName); } + }; + + InitDetInterface(CbmMvdTrackingInterface::Instance(), L1DetectorID::kMvd); + InitDetInterface(CbmStsTrackingInterface::Instance(), L1DetectorID::kSts); + InitDetInterface(CbmMuchTrackingInterface::Instance(), L1DetectorID::kMuch); + InitDetInterface(CbmTrdTrackingInterface::Instance(), L1DetectorID::kTrd); + InitDetInterface(CbmTofTrackingInterface::Instance(), L1DetectorID::kTof); + + InitHitBranch("MvdHit", L1DetectorID::kMvd); + InitHitBranch("StsHit", L1DetectorID::kSts); + InitHitBranch("MuchPixelHit", L1DetectorID::kMuch); + InitHitBranch("TrdHit", L1DetectorID::kTrd); + InitHitBranch("TofHit", L1DetectorID::kTof); + + InitPointBranch("MvdPoint", L1DetectorID::kMvd); + InitPointBranch("StsPoint", L1DetectorID::kSts); + InitPointBranch("MuchPoint", L1DetectorID::kMuch); + InitPointBranch("TrdPoint", L1DetectorID::kTrd); + InitPointBranch("TofPoint", L1DetectorID::kTof); + + this->CheckInit(); + + LOG(error) << fName << ": initializing... \033[1;32mDone\033[0m"; + return kSUCCESS; +} +catch (const std::logic_error& error) { + LOG(error) << fName << ": initializing... \033[1;31mFailed\033[0m\nReason: " << error.what(); + return kFATAL; +} + +// --------------------------------------------------------------------------------------------------------------------- +// +InitStatus InputQaSetup::InitHistograms() +{ + int nStGeo = fpParameters->GetNstationsGeometry(); + + MakeQaDirectory("hit_occupancy"); + for (int iD = 0; iD < static_cast<int>(L1DetectorID::kEND); ++iD) { + if (fvbUseDet[iD]) { MakeQaDirectory(Form("hit_occupancy/%s", kDetName[iD])); } + } + + fph_hit_xz.resize(nStGeo + 1); + fph_hit_yz.resize(nStGeo + 1); + + /// TEMPORARY + constexpr int kNbinsZ = 300; + constexpr int kNbinsX = 200; + constexpr int kNbinsY = 200; + constexpr float kMinZ = -5.; + constexpr float kMaxZ = 350.; + constexpr float kMinX = -100.; + constexpr float kMaxX = 100.; + constexpr float kMinY = -100.; + constexpr float kMaxY = 100.; + + for (int iStGeo = 0; iStGeo <= nStGeo; ++iStGeo) { + auto [detID, iStLoc] = fpParameters->GetStationIndexLocal(iStGeo); + TString sF = ""; + TString sN = ""; + TString sT = ""; + TString sNsuff = (iStGeo == nStGeo) ? "" : Form("_st_%s%d", kDetName[detID], iStLoc); + TString sTsuff = (iStGeo == nStGeo) ? "" : Form(" in %s station %d", kDetName[detID], iStLoc); + + // Hit occupancy + sF = (iStGeo == nStGeo) ? "hit_occupancy" : TString(Form("hit_occupancy/%s", kDetName[detID])); + sN = (TString) "hit_xz" + sNsuff; + sT = (TString) "hit occupancy in xz-plane" + sTsuff + ";z_{hit} [cm];x_{hit} [cm]"; + fph_hit_xz[iStGeo] = MakeQaObject<TH2F>(sF + "/" + sN, sT, kNbinsZ, kMinZ, kMaxZ, kNbinsX, kMinX, kMaxX); + + sN = (TString) "hit_yz" + sNsuff; + sT = (TString) "hit occupancy in yz-plane" + sTsuff + ";z_{hit} [cm];y_{hit} [cm]"; + fph_hit_yz[iStGeo] = MakeQaObject<TH2F>(sF + "/" + sN, sT, kNbinsZ, kMinZ, kMaxZ, kNbinsY, kMinY, kMaxY); + } + + return kSUCCESS; +} + + +// --------------------------------------------------------------------------------------------------------------------- +// +void InputQaSetup::ReadParameters(const char* filename) +{ + fpParameters = std::make_unique<L1Parameters>(); + + L1InitManager manager; + manager.ReadParametersObject(filename); + manager.SendParameters(*fpParameters); + LOG(info) << "CA tracking parameters object: " << filename << '\n' << fpParameters->ToString(1); +} + + +// --------------------------------------------------------------------------------------------------------------------- +// +// --------------------------------------------------------------------------------------------------------------------- +// diff --git a/reco/L1/qa/CbmCaInputQaSetup.h b/reco/L1/qa/CbmCaInputQaSetup.h new file mode 100644 index 0000000000..ee06e67daf --- /dev/null +++ b/reco/L1/qa/CbmCaInputQaSetup.h @@ -0,0 +1,145 @@ +/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt + SPDX-License-Identifier: GPL-3.0-only + Authors: Sergei Zharko [committer] */ + +/// @file CbmCaInputQaSetup.h +/// @brief QA task for tracking detector interfaces +/// @since 28.08.2023 +/// @author S.Zharko <s.zharko@gsi.de> + +#ifndef CbmCaInputQaSetup_h +#define CbmCaInputQaSetup_h 1 + +#include "CbmL1DetectorID.h" +#include "CbmMCDataArray.h" +#include "CbmMCDataObject.h" +#include "CbmMCEventList.h" +#include "CbmMuchPixelHit.h" +#include "CbmMuchPoint.h" +#include "CbmMuchTrackingInterface.h" +#include "CbmMvdHit.h" +#include "CbmMvdPoint.h" +#include "CbmMvdTrackingInterface.h" +#include "CbmQaTask.h" +#include "CbmStsHit.h" +#include "CbmStsPoint.h" +#include "CbmStsTrackingInterface.h" +#include "CbmTimeSlice.h" +#include "CbmTofHit.h" +#include "CbmTofPoint.h" +#include "CbmTofTrackingInterface.h" +#include "CbmTrdHit.h" +#include "CbmTrdPoint.h" +#include "CbmTrdTrackingInterface.h" + +#include "TClonesArray.h" + +#include "L1Parameters.h" + +namespace cbm::ca +{ + /// @class InputQaSetup + /// @brief A QA task to analyze hit and MC point occupancy distributions in different tracking stations + class InputQaSetup : public CbmQaTask { + public: + /// @brief Constructor from parameters + /// @param verbose Verbosity level + /// @param isMCUsed Flag, if MC information is available for this task + InputQaSetup(int verbose, bool isMCUsed); + + /// @brief Reads defined parameters object from file + /// @param filename Name of parameter file + /// @note TEMPORARY FUNCTION, A SEPARATE PARAMETERS INITIALIZATION CLASS IS TO BE USED + void ReadParameters(const char* filename); + + /// @brief Sets detector flag + void SetDetectorFlag(L1DetectorID detID, bool flag = true) { fvbUseDet[detID] = flag; } + + protected: + /// @brief Checks results of the QA and returns a success flag + bool Check(); + + /// @brief Initializes canvases + InitStatus InitCanvases() { return kSUCCESS; }; + + /// @brief Initializes data branches + InitStatus InitDataBranches(); + + /// @brief Initializes histograms + InitStatus InitHistograms(); + + /// @brief Initializes time slice + InitStatus InitTimeSlice() { return kSUCCESS; } + + /// @brief Fills histograms + void FillHistograms(); + + private: + /// @brief Checks branches initialization + void CheckInit() const; + + /// @brief Fill histograms for a given detector type + template<L1DetectorID DetID> + void FillHistogramsDet(); + + + // ******************************* + // ** Data branches and flags ** + // ******************************* + /// @brief Pointers to the tracking detector interfaces for each subsystem + DetIdArr_t<const CbmTrackingDetectorInterfaceBase*> fvpDetInterface = {{nullptr}}; + + DetIdArr_t<bool> fvbUseDet = {{false}}; ///< Detector subsystem usage flag + DetIdArr_t<TClonesArray*> fvpBrHits = {{nullptr}}; ///< Input branch for hits + DetIdArr_t<CbmMCDataArray*> fvpBrPoints = {{nullptr}}; ///< Input branch for MC points + + CbmMCEventList* fpMCEventList = nullptr; ///< Pointer to MC event list + CbmMCDataObject* fpMCEventHeader = nullptr; ///< Pointer to MC event header + std::unique_ptr<L1Parameters> fpParameters = nullptr; ///< Pointer to CA parameters object + + + // ****************** + // ** Histograms ** + // ****************** + + std::vector<TH2F*> fph_hit_xz = {}; ///< Hit occupancy in x-z plane vs. station + std::vector<TH2F*> fph_hit_yz = {}; ///< Hit occupancy in y-z plane vs. station + std::vector<TH2F*> fph_point_xz = {}; ///< MC point occupancy in x-z plane vs. station + std::vector<TH2F*> fph_point_yz = {}; ///< MC point occupancy in y-z plane vs. station + }; +} // namespace cbm::ca + +// ------------------------------------------------------------------------------------------------------------------- +// +template<L1DetectorID DetID> +void cbm::ca::InputQaSetup::FillHistogramsDet() +{ + using Hit_t = HitTypes_t::at<DetID>; + int nHits = fvpBrHits[DetID]->GetEntriesFast(); + + for (int iH = 0; iH < nHits; ++iH) { + const auto* pHit = dynamic_cast<const Hit_t*>(fvpBrHits[DetID]->At(iH)); + if (!pHit) { + LOG(warn) << fName << ": hit with iH = " << iH << " for detector " << kDetName[DetID] << " is not found"; + } + auto address = pHit->GetAddress(); + + // skip T0 hits + if constexpr (L1DetectorID::kTof == DetID) { + if (5 == CbmTofAddress::GetSmType(address)) { continue; } + } + + int iStLoc = fvpDetInterface[DetID]->GetTrackingStationIndex(address); + int iStGeo = fpParameters->GetStationIndexGeometry(iStLoc, DetID); + auto xHit = pHit->GetX(); + auto yHit = pHit->GetY(); + auto zHit = pHit->GetZ(); + + fph_hit_xz[iStGeo]->Fill(zHit, xHit); + fph_hit_yz[iStGeo]->Fill(zHit, yHit); + fph_hit_xz.back()->Fill(zHit, xHit); + fph_hit_yz.back()->Fill(zHit, yHit); + } // iH +} + +#endif // CbmCaInputQaSetup_h -- GitLab