CA: QA task for setup in tracking

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: []

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']

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

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());

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

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
   ///

reco/L1/qa/CbmCaInputQaSetup.cxx

+/* 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);
+}
+
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+// ---------------------------------------------------------------------------------------------------------------------
+//

reco/L1/qa/CbmCaInputQaSetup.h

+/* 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