From dc62953a50000f2f62fd2acb3546d2ff93aebb5c Mon Sep 17 00:00:00 2001
From: "s.zharko@gsi.de" <s.zharko@gsi.de>
Date: Thu, 6 Apr 2023 13:29:36 +0200
Subject: [PATCH] CA: Output QA updates
- track type and track fit histogrammers added
- index bug in TimeSliceReader class fixed
---
core/qa/CMakeLists.txt | 6 +-
core/qa/CbmQaConstants.h | 23 ++++
core/qa/CbmQaIO.cxx | 36 +++++
core/qa/CbmQaIO.h | 218 +++++++++++++++++++++++++++++++
core/qa/CbmQaTask.cxx | 10 +-
core/qa/CbmQaTask.h | 137 +------------------
macro/run/qa_config.cbm.yaml | 22 ++++
macro/run/run_qa.C | 4 +
reco/L1/CMakeLists.txt | 3 +
reco/L1/CbmCaMCModule.cxx | 15 ++-
reco/L1/CbmCaTimeSliceReader.cxx | 170 +++++++++++++-----------
reco/L1/CbmCaTimeSliceReader.h | 21 +--
reco/L1/CbmL1Hit.h | 19 +++
reco/L1/CbmL1MCTrack.cxx | 67 ++++++++++
reco/L1/CbmL1MCTrack.h | 4 +
reco/L1/CbmL1Track.cxx | 60 +++++++++
reco/L1/CbmL1Track.h | 32 ++++-
reco/L1/L1Algo/L1Constants.h | 22 +++-
reco/L1/L1Algo/L1Field.cxx | 4 +-
reco/L1/L1Algo/L1Fit.h | 2 +-
reco/L1/L1Algo/L1TrackPar.h | 62 ++++++++-
reco/L1/catools/CaToolsMCPoint.h | 203 +++++++++++++++++-----------
reco/L1/qa/CbmCaOutputQa.cxx | 129 ++++++++++++++----
reco/L1/qa/CbmCaOutputQa.h | 67 +++++++++-
reco/L1/qa/CbmCaTrackFitQa.cxx | 122 +++++++++++++++++
reco/L1/qa/CbmCaTrackFitQa.h | 163 +++++++++++++++++++++++
reco/L1/qa/CbmCaTrackTypeQa.cxx | 157 ++++++++++++++++++++++
reco/L1/qa/CbmCaTrackTypeQa.h | 209 +++++++++++++++++++++++++++++
28 files changed, 1637 insertions(+), 350 deletions(-)
create mode 100644 core/qa/CbmQaConstants.h
create mode 100644 core/qa/CbmQaIO.cxx
create mode 100644 core/qa/CbmQaIO.h
create mode 100644 reco/L1/CbmL1Track.cxx
create mode 100644 reco/L1/qa/CbmCaTrackFitQa.cxx
create mode 100644 reco/L1/qa/CbmCaTrackFitQa.h
create mode 100644 reco/L1/qa/CbmCaTrackTypeQa.cxx
create mode 100644 reco/L1/qa/CbmCaTrackTypeQa.h
diff --git a/core/qa/CMakeLists.txt b/core/qa/CMakeLists.txt
index f71105c8a6..d4a6f8e18c 100644
--- a/core/qa/CMakeLists.txt
+++ b/core/qa/CMakeLists.txt
@@ -9,8 +9,12 @@ set(SRCS
CbmQaHist.cxx
CbmQaTable.cxx
CbmQaTask.cxx
+ CbmQaIO.cxx
)
+set(HEADERS
+ CbmQaConstants.h
+ )
set(LIBRARY_NAME CbmQaBase)
set(LINKDEF ${LIBRARY_NAME}LinkDef.h)
@@ -29,6 +33,6 @@ set(INTERFACE_DEPENDENCIES
generate_cbm_library()
-Install(FILES CbmQaTask.h CbmQaCanvas.h CbmQaTable.h CbmQaHist.h CbmQaEff.h CbmQaPie.h
+Install(FILES CbmQaTask.h CbmQaCanvas.h CbmQaTable.h CbmQaHist.h CbmQaEff.h CbmQaPie.h CbmQaConstants.h
DESTINATION include
)
diff --git a/core/qa/CbmQaConstants.h b/core/qa/CbmQaConstants.h
new file mode 100644
index 0000000000..91f62f272f
--- /dev/null
+++ b/core/qa/CbmQaConstants.h
@@ -0,0 +1,23 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// @file CbmQaConstants.h
+/// @brief List of constant expressions, common in QA task
+/// @since 04.04.2023
+/// @author Sergei Zharko <s.zharko@gsi.de>
+
+namespace CbmQaConstants
+{
+
+ // Unit
+ namespace unit
+ {
+ constexpr double cm = 1.;
+ constexpr double m = cm * 1.e+2;
+ constexpr double mm = cm * 1.e-1;
+ constexpr double um = cm * 1.e-4;
+ constexpr double nm = cm * 1.e-7;
+ } // namespace unit
+
+} // namespace CbmQaConstants
\ No newline at end of file
diff --git a/core/qa/CbmQaIO.cxx b/core/qa/CbmQaIO.cxx
new file mode 100644
index 0000000000..da8aac99a5
--- /dev/null
+++ b/core/qa/CbmQaIO.cxx
@@ -0,0 +1,36 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// @file CbmQaIO.cxx
+/// @brief Module for ROOT objects IO interface (implementation)
+/// @author S.Zharko <s.zharko@gsi.de>
+/// @since 29.03.2023
+
+#include "CbmQaIO.h"
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+CbmQaIO::CbmQaIO(const char* folderName, const char* prefixName, TFolder* pParentFolder)
+ : fsRootFolderName(folderName)
+ , fsPrefix(prefixName)
+ , fpParentFolder(pParentFolder)
+{
+ if (fpParentFolder) { fpFolderRoot = fpParentFolder->AddFolder(fsRootFolderName, fsRootFolderName); }
+ else {
+ fpFolderRoot =
+ new TFolder(fsRootFolderName, fsRootFolderName); // The name of the folder follows the name of the task
+ fpFolderRoot->SetOwner(true); // When true, TFolder owns all added objects
+ }
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+CbmQaIO::~CbmQaIO()
+{
+ // Free memory for fpFolderRoot
+ if (fpFolderRoot && !fpParentFolder) {
+ delete fpFolderRoot;
+ fpFolderRoot = nullptr;
+ }
+}
diff --git a/core/qa/CbmQaIO.h b/core/qa/CbmQaIO.h
new file mode 100644
index 0000000000..b6ea8acdb5
--- /dev/null
+++ b/core/qa/CbmQaIO.h
@@ -0,0 +1,218 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// @file CbmQaIO.h
+/// @brief Module for ROOT objects IO interface (header)
+/// @author S.Zharko <s.zharko@gsi.de>
+/// @since 29.03.2023
+
+#ifndef CbmQaIO_h
+#define CbmQaIO_h 1
+
+#include "CbmQaTable.h"
+
+#include "Logger.h"
+
+#include "TCanvas.h"
+#include "TEfficiency.h"
+#include "TFolder.h"
+#include "TH1.h"
+#include "TH2.h"
+#include "TH3.h"
+#include "TParameter.h"
+#include "TProfile.h"
+#include "TProfile2D.h"
+#include "TProfile3D.h"
+#include "TROOT.h"
+#include "TString.h"
+
+/// @brief ROOT object IO interface for QA
+///
+/// The class provides interface to write ROOT object into resulted files
+class CbmQaIO {
+public:
+ enum class EStoringMode
+ {
+ kSAMEDIR, ///< Objects will be stored to root directory
+ kSUBDIR ///< Objects will be stored to corresponding subdirectory
+ };
+
+ /// @brief Constructor
+ /// @param folderName Name of the root folder
+ /// @param prefixName Name of the unique prefix
+ /// @param pParentFolder Pointer to parent folder
+ CbmQaIO(const char* folderName, const char* prefixName, TFolder* pParentFolder = nullptr);
+
+ /// @brief Destructor
+ virtual ~CbmQaIO();
+
+ /// @brief Copy constructor
+ CbmQaIO(const CbmQaIO&) = delete;
+
+ /// @brief Move constructor
+ CbmQaIO(CbmQaIO&&) = delete;
+
+ /// @brief Copy assignment operator
+ CbmQaIO& operator=(const CbmQaIO&) = delete;
+
+ /// @brief Move assignment operator
+ CbmQaIO& operator=(CbmQaIO&&) = delete;
+
+ /// @brief Creates, initializes and registers a canvas
+ /// @tparam T Type of the canvas: TCanvas or CbmQaCanvas (\note T should not be a pointer)
+ /// @param name Name of canvas
+ /// @param args The rest of the arguments, which will be passed to the canvas constructor
+ template<typename T, typename... Args>
+ T* MakeCanvas(const char* name, Args... args);
+
+ /// @brief Creates, initializes and registers an efficiency
+ /// @tparam T Type of the canvas: either TProfile, TProfile1D or TEfficiency (\note T should not be a pointer)
+ /// @param name Name of canvas
+ /// @param args The rest of the arguments, which will be passed to the efficiency constructor
+ template<typename T, typename... Args>
+ T* MakeEfficiency(const char* name, Args... args);
+
+ /// @brief Creates, initializes and registers a histogram
+ /// @tparam T Type of the histogram (\note T should not be a pointer)
+ /// @param name Name of histogram
+ /// @param args The rest of the arguments, which will be passed to the histogram constructor
+ template<typename T, typename... Args>
+ T* MakeHisto(const char* name, Args... args);
+
+ /// @brief Creates, initializes and registers a table
+ /// @param name Name of the table
+ /// @param args The rest of the arguments, which will be passed to the table constructor
+ template<typename... Args>
+ CbmQaTable* MakeTable(const char* name, Args... args);
+
+protected:
+ TString fsRootFolderName = ""; ///< Name of root folder
+ TString fsPrefix = ""; ///< Unique prefix for all writeable root
+
+ EStoringMode fStoringMode = EStoringMode::kSUBDIR; ///< Objects storing mode
+ TFolder* fpFolderRoot = nullptr; ///< Root folder to store histograms and canvases
+ TFolder* fpFolderHist = nullptr; ///< Folder for raw histograms
+ TFolder* fpFolderCanv = nullptr; ///< Folder for canvases
+ TFolder* fpFolderEff = nullptr; ///< Folder for efficiencies
+ TFolder* fpFolderTable = nullptr; ///< Folder for tables
+ TFolder* fpParentFolder = nullptr; ///< Pointer to parent folder
+};
+
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+template<typename T, typename... Args>
+T* CbmQaIO::MakeCanvas(const char* nameBase, Args... args)
+{
+ TString name = fsPrefix + "_" + nameBase;
+ if (gROOT->FindObject(name)) {
+ LOG(warn) << fsRootFolderName << ": A previously created canvas \"" << name << "\" will be deleted";
+ T* pCanv = (T*) gROOT->FindObject(name);
+ delete pCanv;
+ }
+
+ // Create a new canvas
+ T* pCanv = new T(name, args...);
+ pCanv->SetLeftMargin(-1.2);
+ pCanv->SetBottomMargin(-1.2);
+
+
+ // Register canvas in the folder
+ if (fStoringMode == EStoringMode::kSUBDIR) {
+ if (!fpFolderCanv) { fpFolderCanv = fpFolderRoot->AddFolder("canvases", "Canvases"); }
+ fpFolderCanv->Add(pCanv);
+ }
+ else if (fStoringMode == EStoringMode::kSAMEDIR) {
+ fpFolderRoot->Add(pCanv);
+ }
+
+ return pCanv;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+template<typename T, typename... Args>
+T* CbmQaIO::MakeEfficiency(const char* nameBase, Args... args)
+{
+ TString name = fsPrefix + "_" + nameBase;
+ if (gROOT->FindObject(name)) {
+ LOG(warn) << fsRootFolderName << ": A previously created histogram \"" << name << "\" will be deleted";
+ auto* pEff = (T*) gROOT->FindObject(name);
+ delete pEff;
+ }
+
+ // Create a new efficiency
+ auto* pEff = new T(name, args...);
+
+ // Register efficiency in the folder
+ if (fStoringMode == EStoringMode::kSUBDIR) {
+ if (!fpFolderEff) { fpFolderEff = fpFolderRoot->AddFolder("efficiencies", "Efficiencies"); }
+ fpFolderEff->Add(pEff);
+ }
+ else if (fStoringMode == EStoringMode::kSAMEDIR) {
+ fpFolderRoot->Add(pEff);
+ }
+
+ return pEff;
+}
+
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+template<typename T, typename... Args>
+T* CbmQaIO::MakeHisto(const char* nameBase, Args... args)
+{
+ TString name = fsPrefix + "_" + nameBase;
+ // Check, if the histogram with a given name was already created. If so, delete it
+ if (gROOT->FindObject(name)) {
+ LOG(warn) << fsRootFolderName << ": A previously created histogram \"" << name << "\" will be deleted";
+ T* pHist = (T*) gROOT->FindObject(name);
+ delete pHist;
+ }
+
+ T* pHist = new T(name, args...);
+ pHist->SetStats(kTRUE);
+ pHist->Sumw2();
+
+ // Register histogram in the folder
+ if (fStoringMode == EStoringMode::kSUBDIR) {
+ if (!fpFolderHist) { fpFolderHist = fpFolderRoot->AddFolder("histograms", "Histograms"); }
+ fpFolderHist->Add(pHist);
+ }
+ else if (fStoringMode == EStoringMode::kSAMEDIR) {
+ fpFolderRoot->Add(pHist);
+ }
+
+ return pHist;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+template<typename... Args>
+CbmQaTable* CbmQaIO::MakeTable(const char* nameBase, Args... args)
+{
+ TString name = fsPrefix + "_" + nameBase;
+ // Check, if the table with a given name was already created. If so, delete it
+ if (gROOT->FindObject(name)) {
+ LOG(warn) << fsRootFolderName << ": A previously created table \"" << name << "\" will be deleted";
+ CbmQaTable* pTable = (CbmQaTable*) gROOT->FindObject(name);
+ delete pTable;
+ }
+
+ // Create a new table
+ CbmQaTable* pTable = new CbmQaTable(name, args...);
+
+ // Register table in folder
+ if (fStoringMode == EStoringMode::kSUBDIR) {
+ if (!fpFolderTable) { fpFolderTable = fpFolderRoot->AddFolder("tables", "Tables"); }
+ fpFolderTable->Add(pTable);
+ }
+ else if (fStoringMode == EStoringMode::kSAMEDIR) {
+ fpFolderRoot->Add(pTable);
+ }
+
+ return pTable;
+}
+
+#endif // CbmQaIO_h
diff --git a/core/qa/CbmQaTask.cxx b/core/qa/CbmQaTask.cxx
index aed9e84f6c..d01349387c 100644
--- a/core/qa/CbmQaTask.cxx
+++ b/core/qa/CbmQaTask.cxx
@@ -26,9 +26,10 @@ ClassImp(CbmQaTask);
//
CbmQaTask::CbmQaTask(const char* name, const char* prefix, int verbose, bool isMCUsed)
: FairTask(name, verbose)
- , fsPrefix(prefix)
+ , CbmQaIO(name, prefix)
, fbUseMC(isMCUsed)
{
+ fStoringMode = CbmQaIO::EStoringMode::kSUBDIR; // mode of objects arrangement in the output file
}
// ---------------------------------------------------------------------------------------------------------------------
@@ -57,7 +58,7 @@ void CbmQaTask::Finish()
// Write the root folder to sinker
FairSink* pSink = FairRootManager::Instance()->GetSink();
LOG_IF(fatal, !pSink) << fName << ": sink file was not found";
- pSink->WriteObject(fpFolderRoot.get(), nullptr);
+ pSink->WriteObject(fpFolderRoot, nullptr);
}
// ---------------------------------------------------------------------------------------------------------------------
@@ -91,8 +92,6 @@ InitStatus CbmQaTask::Init()
res = std::max(res, InitDataBranches());
// ----- Initialize I/O
- fpFolderRoot = std::make_unique<TFolder>(fName, fName); // The name of the folder follows the name of the task
- fpFolderRoot->SetOwner(true); // When true, TFolder owns all added objects
fpFolderRoot->Add(&fNofEvents);
// ----- Initialize histograms
@@ -121,9 +120,6 @@ void CbmQaTask::DeInitBase()
// De-initialize basic data members
// ...
- delete fpFolderRoot.get();
- fpFolderRoot = nullptr;
-
// De-initialize particular QA-task implementation
DeInit();
}
diff --git a/core/qa/CbmQaTask.h b/core/qa/CbmQaTask.h
index 143b7a7abd..7f9e39e71b 100644
--- a/core/qa/CbmQaTask.h
+++ b/core/qa/CbmQaTask.h
@@ -5,11 +5,12 @@
/// \file CbmQaTask.h
/// \brief A base class for CBM QA task logic
/// \author S.Zharko <s.zharko@gsi.de>
-/// \data 12.01.2022
+/// \data 12.01.2023
#ifndef CbmQaTask_h
#define CbmQaTask_h 1
+#include "CbmQaIO.h"
#include "CbmQaTable.h"
#include "FairTask.h"
@@ -36,7 +37,7 @@
/// Class CbmQaTask is to be inherited with a particular QA-task. It provides mechanisms for storage and management
/// of QA canvases and histograms management
-class CbmQaTask : public FairTask {
+class CbmQaTask : public FairTask, public CbmQaIO {
public:
/// Constructor from parameters
/// \param name Name of the task
@@ -108,27 +109,6 @@ protected:
/// Method to fill histograms per event or time-slice
virtual void FillHistograms() {}
- /// Creates, initializes and registers a canvas
- /// \tparam T Type of the canvas: TCanvas or CbmQaCanvas (\note T should not be a pointer)
- /// \param name Name of canvas
- /// \param args The rest of the arguments, which will be passed to the canvas constructor
- template<typename T, typename... Args>
- T* MakeCanvas(const char* name, Args... args);
-
- /// Creates, initializes and registers an efficiency
- /// \tparam T Type of the canvas: either TProfile, TProfile2D or TEfficiency (\note T should not be a pointer)
- /// \param name Name of canvas
- /// \param args The rest of the arguments, which will be passed to the efficiency constructor
- template<typename T, typename... Args>
- T* MakeEfficiency(const char* name, Args... args);
-
- /// Creates, initializes and registers a histogram
- /// \tparam T Type of the histogram (\note T should not be a pointer)
- /// \param name Name of histogram
- /// \param args The rest of the arguments, which will be passed to the histogram constructor
- template<typename T, typename... Args>
- T* MakeHisto(const char* name, Args... args);
-
/// @brief Creates, initializes and registers a histogram, based on the configuration file
/// @tparam T Type of the histogram (@note: should not be a pointer)
/// @param name Name of the histogram, stored in config
@@ -139,12 +119,6 @@ protected:
template<typename T>
T* MakeHistoFromConfig(const char* name, int id0 = -1, int id1 = -1, int id2 = -1);
- /// Creates, initializes and registers a table
- /// \param name Name of the table
- /// \param args The rest of the arguments, which will be passed to the table constructor
- template<typename... Args>
- CbmQaTable* MakeTable(const char* name, Args... args);
-
/// Get current event number
int GetEventNumber() const { return fNofEvents.GetVal(); }
@@ -160,21 +134,14 @@ protected:
/// \param hi Upper limit of the variable
/// \return False, if variable exits the range
template<typename T>
- bool CheckRange(const std::string_view& name, const T& var, const T& lo, const T& hi) const;
+ bool CheckRange(std::string_view name, const T& var, const T& lo, const T& hi) const;
private:
/// @brief De-initializes this task
void DeInitBase();
- // I/O
- std::unique_ptr<TFolder> fpFolderRoot = nullptr; ///< Root folder to store histograms and canvases
- TFolder* fpFolderHist = nullptr; ///< Folder for raw histograms
- TFolder* fpFolderCanv = nullptr; ///< Folder for canvases
- TFolder* fpFolderEff = nullptr; ///< Folder for efficiencies
- TFolder* fpFolderTable = nullptr; ///< Folder for tables
- TString fsPrefix = ""; ///< Unique prefix for all writeable root objects to avoid collisions between different tasks
bool fbUseMC = false; ///< Flag, if MC is used
TParameter<int> fNofEvents {"nEvents", 0}; ///< Number of processed events
@@ -191,7 +158,7 @@ private:
// ---------------------------------------------------------------------------------------------------------------------
//
template<typename T>
-bool CbmQaTask::CheckRange(const std::string_view& name, const T& var, const T& lo, const T& hi) const
+bool CbmQaTask::CheckRange(std::string_view name, const T& var, const T& lo, const T& hi) const
{
if (std::clamp(var, lo, hi) == lo) {
LOG(error) << fName << ": " << name << " is found to be under the lower limit (" << var << " < " << lo << ')';
@@ -204,78 +171,6 @@ bool CbmQaTask::CheckRange(const std::string_view& name, const T& var, const T&
return true;
}
-// ---------------------------------------------------------------------------------------------------------------------
-//
-template<typename T, typename... Args>
-T* CbmQaTask::MakeCanvas(const char* nameBase, Args... args)
-{
- TString name = fsPrefix + "_" + nameBase;
- if (gROOT->FindObject(name)) {
- LOG(warn) << fName << ": A canvas with name \"" << name << "\" was previously created and will be deleted now "
- << "to avoid memory leaks";
- T* pCanv = (T*) gROOT->FindObject(name);
- delete pCanv;
- }
-
- // Create a new canvas
- T* pCanv = new T(name, args...);
- pCanv->SetLeftMargin(0.2);
- pCanv->SetBottomMargin(0.2);
-
-
- // Register canvas in the folder
- if (!fpFolderCanv) { fpFolderCanv = fpFolderRoot->AddFolder("canvases", "Canvases"); }
- fpFolderCanv->Add(pCanv);
-
- return pCanv;
-}
-
-// ---------------------------------------------------------------------------------------------------------------------
-//
-template<typename T, typename... Args>
-T* CbmQaTask::MakeEfficiency(const char* nameBase, Args... args)
-{
- TString name = fsPrefix + "_" + nameBase;
- if (gROOT->FindObject(name)) {
- LOG(warn) << fName << ": An efficiency with name \"" << name << "\" was previously created and will be deleted now "
- << "to avoid memory leaks";
- auto* pEff = (T*) gROOT->FindObject(name);
- delete pEff;
- }
-
- // Create a new efficiency
- auto* pEff = new T(name, args...);
-
- // Register efficiency in the folder
- if (!fpFolderEff) { fpFolderEff = fpFolderRoot->AddFolder("efficiencies", "Efficiencies"); }
- fpFolderEff->Add(pEff);
-
- return pEff;
-}
-
-
-// ---------------------------------------------------------------------------------------------------------------------
-//
-template<typename T, typename... Args>
-T* CbmQaTask::MakeHisto(const char* nameBase, Args... args)
-{
- TString name = fsPrefix + "_" + nameBase;
- // Check, if the histogram with a given name was already created. If so, delete it
- if (gROOT->FindObject(name)) {
- LOG(warn) << fName << ": A histogram with name \"" << name << "\" was previously created and will be deleted now "
- << "to avoid memory leaks";
- T* pHist = (T*) gROOT->FindObject(name);
- delete pHist;
- }
-
- T* pHist = new T(name, args...);
-
- // Register histogram in the folder
- if (!fpFolderHist) { fpFolderHist = fpFolderRoot->AddFolder("histograms", "Histograms"); }
- fpFolderHist->Add(pHist);
-
- return pHist;
-}
// ---------------------------------------------------------------------------------------------------------------------
//
@@ -397,29 +292,7 @@ T* CbmQaTask::MakeHistoFromConfig(const char* nameBase, int id0, int id1, int id
return pHist;
}
-// ---------------------------------------------------------------------------------------------------------------------
-//
-template<typename... Args>
-CbmQaTable* CbmQaTask::MakeTable(const char* nameBase, Args... args)
-{
- TString name = fsPrefix + "_" + nameBase;
- // Check, if the table with a given name was already created. If so, delete it
- if (gROOT->FindObject(name)) {
- LOG(warn) << fName << ": A table with name \"" << name << "\" was previously created and will be deleted now "
- << "to avoid memory leaks";
- CbmQaTable* pTable = (CbmQaTable*) gROOT->FindObject(name);
- delete pTable;
- }
-
- // Create a new table
- CbmQaTable* pTable = new CbmQaTable(name, args...);
-
- // Register table in folder
- if (!fpFolderTable) { fpFolderTable = fpFolderRoot->AddFolder("tables", "Tables"); }
- fpFolderTable->Add(pTable);
- return pTable;
-}
#endif // CbmQaTask_h
diff --git a/macro/run/qa_config.cbm.yaml b/macro/run/qa_config.cbm.yaml
index 6d23fa8722..baf8957e05 100644
--- a/macro/run/qa_config.cbm.yaml
+++ b/macro/run/qa_config.cbm.yaml
@@ -208,6 +208,28 @@ qa:
miny: -0.5
maxy: 14.5
#
+ # -- Reconstructed tracks vs. MC quantities
+ #
+ reco_pMC:
+ title: "MC momentum of reconstructed tracks;p_{MC} [GeV/c]"
+ nbinsx: 150
+ minx: 0.
+ maxx: 15.
+ reco_yMC:
+ title: "MC momentum of reconstructed tracks;p_{MC} [GeV/c]"
+ nbinsy: 40
+ miny: -4.
+ maxy: +4.
+ reco_pMC_vs_yMC:
+ title: "MC momentum vs. MC rapidity of reconstructed tracks;p_{MC} [GeV/c]"
+ nbinsx: 150
+ minx: 0.
+ maxx: 15.
+ nbinsy: 40
+ miny: -4.
+ maxy: +4.
+
+ #
# -- Track residuals and pulls --
#
# Residuals and pulls at first station
diff --git a/macro/run/run_qa.C b/macro/run/run_qa.C
index bb9f108253..296094afc5 100644
--- a/macro/run/run_qa.C
+++ b/macro/run/run_qa.C
@@ -249,6 +249,10 @@ void run_qa(TString dataTra = "data/sis100_muon_jpsi_test", TString dataRaw = "d
// ------------------------------------------------------------------------
// ----- Tracking QA ------------------------------------------------------
+ // KF is currently needed to access magnetic field. In future we will
+ // delegate track fit routines to CbmKF as well.
+ run->AddTask(new CbmKF());
+
TString caParFile = recFile;
caParFile.ReplaceAll(".root", ".L1Parameters.dat");
diff --git a/reco/L1/CMakeLists.txt b/reco/L1/CMakeLists.txt
index ddad3866d3..60ed071a1c 100644
--- a/reco/L1/CMakeLists.txt
+++ b/reco/L1/CMakeLists.txt
@@ -54,6 +54,7 @@ set(SRCS
L1Algo/L1Fit.cxx
L1Algo/L1MCEvent.cxx
CbmL1MCTrack.cxx
+ CbmL1Track.cxx
L1Algo/L1Material.cxx
L1Algo/L1UMeasurementInfo.cxx
L1Algo/L1XYMeasurementInfo.cxx
@@ -88,6 +89,8 @@ set(SRCS
qa/CbmCaInputQaTrd.cxx
qa/CbmCaInputQaTof.cxx
qa/CbmCaOutputQa.cxx
+ qa/CbmCaTrackTypeQa.cxx
+ qa/CbmCaTrackFitQa.cxx
qa/CbmTofInteraction.cxx # Tests
)
diff --git a/reco/L1/CbmCaMCModule.cxx b/reco/L1/CbmCaMCModule.cxx
index 3f56f34819..ba105dc7e5 100644
--- a/reco/L1/CbmCaMCModule.cxx
+++ b/reco/L1/CbmCaMCModule.cxx
@@ -427,7 +427,7 @@ void CbmCaMCModule::MatchRecoAndMCTracks()
// ----- Count number of hits from each MC track belonging to this reconstructed track
auto& mNofHitsVsMCTrkID = trkRe.hitMap;
- mNofHitsVsMCTrkID.clear();
+ //mNofHitsVsMCTrkID.clear();
for (int iH : trkRe.Hits) {
for (int iP : (*fpvQaHits)[iH].mcPointIds) {
assert(iP > -1); // Should never be triggered
@@ -436,7 +436,6 @@ void CbmCaMCModule::MatchRecoAndMCTracks()
else {
mNofHitsVsMCTrkID[iTmc] += 1;
}
-
} // loop over global point ids stored for hit: end
} // loop over hit ids stored for a reconstructed track: end
@@ -470,7 +469,9 @@ void CbmCaMCModule::MatchRecoAndMCTracks()
// Update max purity of the reconstructed track
trkRe.SetMaxPurity(maxPurity);
} // loop over hit map: end
- } // loop over reconstructed tracks: end
+ assert(trkRe.GetNofMCTracks() < 2);
+
+ } // loop over reconstructed tracks: end
}
@@ -695,6 +696,7 @@ void CbmCaMCModule::ReadMCTracks()
if (!pEvtHeader) {
LOG(fatal) << "CbmCaMCModule: event header is not found for file " << iFile << " and event " << iEvent;
}
+ double eventTime = fpMCEventList->GetEventTime(iEvent, iFile);
// ----- Loop over MC tracks
int nTracks = fpMCTracks->Size(iFile, iEvent);
@@ -712,7 +714,7 @@ void CbmCaMCModule::ReadMCTracks()
aTrk.SetEventId(iEvent);
aTrk.SetFileId(iFile);
- aTrk.SetStartT(pExtMCTrk->GetStartT());
+ aTrk.SetStartT(pExtMCTrk->GetStartT() + eventTime);
aTrk.SetStartX(pExtMCTrk->GetStartX());
aTrk.SetStartY(pExtMCTrk->GetStartY());
aTrk.SetStartZ(pExtMCTrk->GetStartZ());
@@ -735,11 +737,12 @@ void CbmCaMCModule::ReadMCTracks()
// Set index of mother track. We assume, that mother track was recorded into the internal array before
int extMotherId = pExtMCTrk->GetMotherId();
if (extMotherId < 0) {
- // This is a primary track, and it's mother ID equals -1 or -2. This value is taken also as internal mother ID
+ // This is a primary track, and it's mother ID equals -1 or -2. This value is taken also as an internal mother
+ // ID.
aTrk.SetMotherId(extMotherId);
}
else {
- // This is a secondary track, mother ID should be recalculated for the internal array
+ // This is a secondary track, mother ID should be recalculated for the internal track array.
int motherId = fpMCData->FindInternalTrackIndex(extMotherId, iEvent, iFile);
if (motherId == -1) { motherId = -3; } // Mother is neutral particle, which is rejected
aTrk.SetMotherId(motherId);
diff --git a/reco/L1/CbmCaTimeSliceReader.cxx b/reco/L1/CbmCaTimeSliceReader.cxx
index 1f55a04df6..63607559f2 100644
--- a/reco/L1/CbmCaTimeSliceReader.cxx
+++ b/reco/L1/CbmCaTimeSliceReader.cxx
@@ -16,6 +16,8 @@
#include "FairRootManager.h"
#include "Logger.h"
+#include <numeric>
+
#include "L1Constants.h"
#include "L1InputData.h"
#include "L1Parameters.h"
@@ -161,76 +163,6 @@ void TimeSliceReader::InitTimeSlice()
this->ReadRecoTracks();
}
-// ---------------------------------------------------------------------------------------------------------------------
-//
-void TimeSliceReader::ReadHits()
-{
- fNofHits = 0;
- fNofHitKeys = 0;
- fFirstHitKey = 0;
-
- // Get total number of hits
- int nHitsTot = 0;
- if (fbUseMvd) { nHitsTot += fpBrMvdHits->GetEntriesFast(); }
- if (fbUseSts) { nHitsTot += fpBrStsHits->GetEntriesFast(); }
- if (fbUseMuch) { nHitsTot += fpBrMuchHits->GetEntriesFast(); }
- if (fbUseTrd) { nHitsTot += fpBrTrdHits->GetEntriesFast(); }
- if (fbUseTof) { nHitsTot += fpBrTofHits->GetEntriesFast(); }
-
- // Resize the containers
- if (fpvHitIds) {
- fpvHitIds->clear();
- fpvHitIds->reserve(nHitsTot);
- }
- if (fpvQaHits) {
- fpvQaHits->clear();
- fpvQaHits->reserve(nHitsTot);
- }
- if (fpIODataManager) {
- fpIODataManager->ResetInputData();
- fpIODataManager->ReserveNhits(nHitsTot);
- }
-
- std::fill(fvHitFirstIndexDet.begin(), fvHitFirstIndexDet.end(), 0);
-
- // Save number of hits stored
- std::array<int, L1Constants::size::kMaxNdetectors> vNofHits = {0};
-
- // Read hits for different detectors
- if (fbUseMvd) { vNofHits[int(L1DetectorID::kMvd)] += ReadHitsForDetector<L1DetectorID::kMvd>(fpBrMvdHits); }
- if (fbUseSts) { vNofHits[int(L1DetectorID::kSts)] += ReadHitsForDetector<L1DetectorID::kSts>(fpBrStsHits); }
- if (fbUseMuch) { vNofHits[int(L1DetectorID::kMuch)] += ReadHitsForDetector<L1DetectorID::kMuch>(fpBrMuchHits); }
- if (fbUseTrd) { vNofHits[int(L1DetectorID::kTrd)] += ReadHitsForDetector<L1DetectorID::kTrd>(fpBrTrdHits); }
- if (fbUseTof) { vNofHits[int(L1DetectorID::kTof)] += ReadHitsForDetector<L1DetectorID::kTof>(fpBrTofHits); }
-
- // Save first hit index for different detector subsystems
- for (uint32_t iDet = 0; iDet < vNofHits.size(); ++iDet) {
- fvHitFirstIndexDet[iDet + 1] = fvHitFirstIndexDet[iDet] + vNofHits[iDet];
- fNofHits += vNofHits[iDet];
- }
-
- // Update maps of ext->int hit indexes
- // NOTE: fvpHitIds must be initialized, if we want to read tracks from the file
- if (fpvHitIds) {
- auto UpdateHitIndexMap = [&](L1Vector<int> & vIds, L1DetectorID detID) constexpr
- {
- vIds.reset(vNofHits[int(detID)]);
- for (int iH = fvHitFirstIndexDet[int(detID)]; iH < fvHitFirstIndexDet[int(detID) + 1]; ++iH) {
- vIds[(*fpvHitIds)[iH].hitId] = iH;
- }
- };
- if (fbUseMvd) { UpdateHitIndexMap(vHitMvdIds, L1DetectorID::kMvd); }
- if (fbUseSts) { UpdateHitIndexMap(vHitStsIds, L1DetectorID::kSts); }
- if (fbUseMuch) { UpdateHitIndexMap(vHitMuchIds, L1DetectorID::kMuch); }
- if (fbUseTrd) { UpdateHitIndexMap(vHitTrdIds, L1DetectorID::kTrd); }
- if (fbUseTof) { UpdateHitIndexMap(vHitTofIds, L1DetectorID::kTof); }
- }
- // Update number of hit keys in input data object
- if (fpIODataManager) { fpIODataManager->SetNhitKeys(fNofHitKeys); }
-
- // Sort debug hits
- if (fpvQaHits) { this->SortQaHits(); }
-}
// ---------------------------------------------------------------------------------------------------------------------
//
@@ -263,11 +195,13 @@ void TimeSliceReader::ReadRecoTracks()
track.Hits.clear();
track.Hits.reserve(pInputTrack->GetNofHits());
for (int iH = 0; iH < pInputTrack->GetNofMvdHits(); ++iH) {
- int iHint = vHitMvdIds[pInputTrack->GetMvdHitIndex(iH)];
- track.Hits.push_back(iHint); // !!!!
- } // iH
+ int iHext = pInputTrack->GetMvdHitIndex(iH);
+ int iHint = vHitMvdIds[iHext];
+ track.Hits.push_back(iHint);
+ } // iH
for (int iH = 0; iH < pInputTrack->GetNofStsHits(); ++iH) {
- int iHint = vHitStsIds[pInputTrack->GetStsHitIndex(iH)];
+ int iHext = pInputTrack->GetStsHitIndex(iH);
+ int iHint = vHitStsIds[iHext];
track.Hits.push_back(iHint);
} // iH
} // iT
@@ -303,7 +237,8 @@ void TimeSliceReader::SetDetector(L1DetectorID detID, bool flag)
// ---------------------------------------------------------------------------------------------------------------------
//
-void TimeSliceReader::SortQaHits()
+template<>
+void TimeSliceReader::SortQaHits<true>()
{
int nStationsActive = fpParameters->GetNstationsActive();
L1Vector<CbmL1HitDebugInfo> vNewHits(fpvQaHits->size());
@@ -328,6 +263,90 @@ void TimeSliceReader::SortQaHits()
std::swap(vNewHits, (*fpvQaHits));
}
+// ---------------------------------------------------------------------------------------------------------------------
+//
+template<>
+void TimeSliceReader::SortQaHits<false>()
+{
+ std::sort(fpvQaHits->begin(), fpvQaHits->end(), [](const CbmL1HitDebugInfo& a, const CbmL1HitDebugInfo& b) {
+ return (a.iStation < b.iStation) || (a.iStation == b.iStation) && (a.y < b.y);
+ });
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void TimeSliceReader::ReadHits()
+{
+ fNofHits = 0;
+ fNofHitKeys = 0;
+ fFirstHitKey = 0;
+
+ // Get total number of hits
+ std::array<int, L1Constants::size::kMaxNdetectors> nHitsDet = {0};
+ if (fbUseMvd) { nHitsDet[int(L1DetectorID::kMvd)] = fpBrMvdHits->GetEntriesFast(); }
+ if (fbUseSts) { nHitsDet[int(L1DetectorID::kSts)] = fpBrStsHits->GetEntriesFast(); }
+ if (fbUseMuch) { nHitsDet[int(L1DetectorID::kMuch)] = fpBrMuchHits->GetEntriesFast(); }
+ if (fbUseTrd) { nHitsDet[int(L1DetectorID::kTrd)] = fpBrTrdHits->GetEntriesFast(); }
+ if (fbUseTof) { nHitsDet[int(L1DetectorID::kTof)] = fpBrTofHits->GetEntriesFast(); }
+
+ int nHitsTot = std::accumulate(nHitsDet.begin(), nHitsDet.end(), 0);
+
+ // Resize the containers
+ if (fpvHitIds) {
+ fpvHitIds->clear();
+ fpvHitIds->reserve(nHitsTot);
+ }
+ if (fpvQaHits) {
+ fpvQaHits->clear();
+ fpvQaHits->reserve(nHitsTot);
+ }
+ if (fpIODataManager) {
+ fpIODataManager->ResetInputData();
+ fpIODataManager->ReserveNhits(nHitsTot);
+ }
+
+ std::fill(fvHitFirstIndexDet.begin(), fvHitFirstIndexDet.end(), 0);
+
+ // Save number of hits stored
+ std::array<int, L1Constants::size::kMaxNdetectors> vNofHits = {0};
+
+ // Read hits for different detectors
+ if (fbUseMvd) { vNofHits[int(L1DetectorID::kMvd)] += ReadHitsForDetector<L1DetectorID::kMvd>(fpBrMvdHits); }
+ if (fbUseSts) { vNofHits[int(L1DetectorID::kSts)] += ReadHitsForDetector<L1DetectorID::kSts>(fpBrStsHits); }
+ if (fbUseMuch) { vNofHits[int(L1DetectorID::kMuch)] += ReadHitsForDetector<L1DetectorID::kMuch>(fpBrMuchHits); }
+ if (fbUseTrd) { vNofHits[int(L1DetectorID::kTrd)] += ReadHitsForDetector<L1DetectorID::kTrd>(fpBrTrdHits); }
+ if (fbUseTof) { vNofHits[int(L1DetectorID::kTof)] += ReadHitsForDetector<L1DetectorID::kTof>(fpBrTofHits); }
+
+ // Save first hit index for different detector subsystems
+ for (uint32_t iDet = 0; iDet < vNofHits.size(); ++iDet) {
+ fvHitFirstIndexDet[iDet + 1] = fvHitFirstIndexDet[iDet] + vNofHits[iDet];
+ fNofHits += vNofHits[iDet];
+ }
+
+ // Update number of hit keys in input data object
+ if (fpIODataManager) { fpIODataManager->SetNhitKeys(fNofHitKeys); }
+
+ // Sort debug hits
+ if (fpvQaHits) { this->SortQaHits<true>(); } // false - old sort, true - new backet sort
+
+ // Update maps of ext->int hit indexes
+ // NOTE: fvpHitIds must be initialized, if we want to read tracks from the file
+ if (fpvHitIds) {
+ auto UpdateHitIndexMap = [&](L1Vector<int> & vIds, L1DetectorID detID) constexpr
+ {
+ vIds.reset(nHitsDet[int(detID)]);
+ for (int iH = fvHitFirstIndexDet[int(detID)]; iH < fvHitFirstIndexDet[int(detID) + 1]; ++iH) {
+ vIds[(*fpvQaHits)[iH].ExtIndex] = iH;
+ }
+ };
+ if (fbUseMvd) { UpdateHitIndexMap(vHitMvdIds, L1DetectorID::kMvd); }
+ if (fbUseSts) { UpdateHitIndexMap(vHitStsIds, L1DetectorID::kSts); }
+ if (fbUseMuch) { UpdateHitIndexMap(vHitMuchIds, L1DetectorID::kMuch); }
+ if (fbUseTrd) { UpdateHitIndexMap(vHitTrdIds, L1DetectorID::kTrd); }
+ if (fbUseTof) { UpdateHitIndexMap(vHitTofIds, L1DetectorID::kTof); }
+ }
+}
+
// ---------------------------------------------------------------------------------------------------------------------
//
void TimeSliceReader::StoreHitRecord(const HitRecord& hitRecord)
@@ -357,6 +376,7 @@ void TimeSliceReader::StoreHitRecord(const HitRecord& hitRecord)
CbmL1HitDebugInfo hitDbg;
hitDbg.Det = hitRecord.fDet;
hitDbg.ExtIndex = hitRecord.fExtId;
+ hitDbg.IntIndex = static_cast<int>(fpvQaHits->size());
hitDbg.iStation = hitRecord.fStaId;
hitDbg.x = hitRecord.fX;
hitDbg.y = hitRecord.fY;
diff --git a/reco/L1/CbmCaTimeSliceReader.h b/reco/L1/CbmCaTimeSliceReader.h
index a92b7809a8..ac4ffa17b7 100644
--- a/reco/L1/CbmCaTimeSliceReader.h
+++ b/reco/L1/CbmCaTimeSliceReader.h
@@ -28,6 +28,7 @@
#include "TClonesArray.h"
+#include "L1Constants.h"
#include "L1Vector.h"
@@ -162,6 +163,7 @@ namespace cbm::ca
int ReadHitsForDetector(const TClonesArray* pBrHits);
/// @brief Sorts QA hit objects by stations
+ template<bool IsNewSortingApproach>
void SortQaHits();
/// @brief Saves hit to data structures
@@ -170,6 +172,7 @@ namespace cbm::ca
/// Stores recorded hit information into registered hit containers
void StoreHitRecord(const HitRecord& hitRecord);
+
// Flags for detector subsystems being used
bool fbUseMvd = false; ///< is MVD used
bool fbUseSts = false; ///< is STS used
@@ -212,6 +215,7 @@ namespace cbm::ca
std::shared_ptr<L1IODataManager> fpIODataManager = nullptr; ///< Pointer to input data manager
std::shared_ptr<L1Parameters> fpParameters = nullptr; ///< Pointer to tracking parameters object
+
// Maps of hit indexes: ext -> int
L1Vector<int> vHitMvdIds {"CbmCaTimeSliceReader::vHitMvdIds"};
L1Vector<int> vHitStsIds {"CbmCaTimeSliceReader::vHitStsIds"};
@@ -219,6 +223,7 @@ namespace cbm::ca
L1Vector<int> vHitTrdIds {"CbmCaTimeSliceReader::vHitTrdIds"};
L1Vector<int> vHitTofIds {"CbmCaTimeSliceReader::vHitTofIds"};
+
// Additional
ECbmTrackingMode fTrackingMode; ///< Tracking mode
@@ -247,7 +252,7 @@ int cbm::ca::TimeSliceReader::ReadHitsForDetector(const TClonesArray* pBrHits)
fFirstHitKey = fNofHitKeys;
- for (int iH = 0; iH < nHitsTot; ++iH) {
+ for (int iHext = 0; iHext < nHitsTot; ++iHext) {
HitRecord hitRecord; // Record of hits information
CbmPixelHit* pPixelHit = nullptr; // Pointer to hit object
@@ -257,7 +262,7 @@ int cbm::ca::TimeSliceReader::ReadHitsForDetector(const TClonesArray* pBrHits)
// Fill out detector specific data
if constexpr (L1DetectorID::kMvd == DetID) {
- CbmMvdHit* pMvdHit = static_cast<CbmMvdHit*>(pBrHits->At(iH));
+ CbmMvdHit* pMvdHit = static_cast<CbmMvdHit*>(pBrHits->At(iHext));
iStGeom = fpMvdInterface->GetTrackingStationIndex(pMvdHit->GetStationNr());
phiF = fpMvdInterface->GetStripsStereoAngleFront(iStGeom);
phiB = fpMvdInterface->GetStripsStereoAngleBack(iStGeom);
@@ -266,7 +271,7 @@ int cbm::ca::TimeSliceReader::ReadHitsForDetector(const TClonesArray* pBrHits)
hitRecord.fDv = pMvdHit->GetDy();
}
else if constexpr (L1DetectorID::kSts == DetID) {
- CbmStsHit* pStsHit = static_cast<CbmStsHit*>(pBrHits->At(iH));
+ CbmStsHit* pStsHit = static_cast<CbmStsHit*>(pBrHits->At(iHext));
iStGeom = fpStsInterface->GetTrackingStationIndex(pStsHit->GetAddress());
phiF = fpStsInterface->GetStripsStereoAngleFront(iStGeom);
phiB = fpStsInterface->GetStripsStereoAngleBack(iStGeom);
@@ -277,7 +282,7 @@ int cbm::ca::TimeSliceReader::ReadHitsForDetector(const TClonesArray* pBrHits)
hitRecord.fDv = pStsHit->GetDv();
}
else if constexpr (L1DetectorID::kMuch == DetID) {
- CbmMuchPixelHit* pMuchHit = static_cast<CbmMuchPixelHit*>(pBrHits->At(iH));
+ CbmMuchPixelHit* pMuchHit = static_cast<CbmMuchPixelHit*>(pBrHits->At(iHext));
iStGeom = fpMuchInterface->GetTrackingStationIndex(pMuchHit->GetAddress());
phiF = fpMuchInterface->GetStripsStereoAngleFront(iStGeom);
phiB = fpMuchInterface->GetStripsStereoAngleBack(iStGeom);
@@ -286,7 +291,7 @@ int cbm::ca::TimeSliceReader::ReadHitsForDetector(const TClonesArray* pBrHits)
hitRecord.fDv = pMuchHit->GetDy();
}
else if constexpr (L1DetectorID::kTrd == DetID) {
- CbmTrdHit* pTrdHit = static_cast<CbmTrdHit*>(pBrHits->At(iH));
+ CbmTrdHit* pTrdHit = static_cast<CbmTrdHit*>(pBrHits->At(iHext));
iStGeom = fpTrdInterface->GetTrackingStationIndex(pTrdHit->GetAddress());
phiF = fpTrdInterface->GetStripsStereoAngleFront(iStGeom);
phiB = fpTrdInterface->GetStripsStereoAngleBack(iStGeom);
@@ -295,7 +300,7 @@ int cbm::ca::TimeSliceReader::ReadHitsForDetector(const TClonesArray* pBrHits)
hitRecord.fDv = pTrdHit->GetDy();
}
else if constexpr (L1DetectorID::kTof == DetID) {
- CbmTofHit* pTofHit = static_cast<CbmTofHit*>(pBrHits->At(iH));
+ CbmTofHit* pTofHit = static_cast<CbmTofHit*>(pBrHits->At(iHext));
iStGeom = fpTofInterface->GetTrackingStationIndex(pTofHit->GetAddress());
phiF = fpTofInterface->GetStripsStereoAngleFront(iStGeom);
phiB = fpTofInterface->GetStripsStereoAngleBack(iStGeom);
@@ -328,7 +333,7 @@ int cbm::ca::TimeSliceReader::ReadHitsForDetector(const TClonesArray* pBrHits)
hitRecord.fDataStream = (static_cast<int64_t>(hitRecord.fDet) << 60) | pPixelHit->GetAddress();
hitRecord.fU = hitRecord.fX * cos(phiF) + hitRecord.fY * sin(phiF);
hitRecord.fV = hitRecord.fX * cos(phiB) + hitRecord.fY * sin(phiB);
- hitRecord.fExtId = iH;
+ hitRecord.fExtId = iHext;
// Update number of hit keys
if constexpr (L1DetectorID::kSts == DetID) {
@@ -336,7 +341,7 @@ int cbm::ca::TimeSliceReader::ReadHitsForDetector(const TClonesArray* pBrHits)
if (fNofHitKeys <= hitRecord.fStripB) { fNofHitKeys += hitRecord.fStripB; }
}
else {
- hitRecord.fStripF = fFirstHitKey + iH;
+ hitRecord.fStripF = fFirstHitKey + iHext;
hitRecord.fStripB = hitRecord.fStripF;
if (fNofHitKeys <= hitRecord.fStripF) { fNofHitKeys += hitRecord.fStripF; }
}
diff --git a/reco/L1/CbmL1Hit.h b/reco/L1/CbmL1Hit.h
index 462548ffce..e6a2901336 100644
--- a/reco/L1/CbmL1Hit.h
+++ b/reco/L1/CbmL1Hit.h
@@ -23,6 +23,22 @@ public:
CbmL1HitId() = default;
CbmL1HitId(int det, int index) : detId(det), hitId(index) {};
+ /// @brief String representation of class object
+ /// @param header If true, header will be printed
+ std::string ToString(bool header = false) const
+ {
+ std::stringstream msg;
+ if (header) {
+ msg << std::setw(8) << std::setfill(' ') << "det. id" << ' ';
+ msg << std::setw(8) << std::setfill(' ') << "ext. id";
+ }
+ else {
+ msg << std::setw(8) << std::setfill(' ') << detId << ' ';
+ msg << std::setw(8) << std::setfill(' ') << hitId;
+ }
+ return msg.str();
+ }
+
int detId = -1; ///< detector ID
int hitId = -1; ///< index of hit in the TClonesArray array
};
@@ -76,6 +92,7 @@ public:
std::stringstream msg;
if (header) {
msg << setw(8) << setfill(' ') << "ext. ID" << ' ';
+ msg << setw(8) << setfill(' ') << "int. ID" << ' ';
msg << setw(8) << setfill(' ') << "st. ID" << ' ';
msg << setw(12) << setfill(' ') << "x [cm]" << ' ';
msg << setw(12) << setfill(' ') << "y [cm]" << ' ';
@@ -84,6 +101,7 @@ public:
}
else {
msg << setw(8) << setfill(' ') << ExtIndex << ' ';
+ msg << setw(8) << setfill(' ') << IntIndex << ' ';
msg << setw(8) << setfill(' ') << iStation << ' ';
msg << setw(12) << setfill(' ') << x << ' ';
msg << setw(12) << setfill(' ') << y << ' ';
@@ -101,6 +119,7 @@ public:
// TODO: SZh 2.03.2023: make the variables private
int ExtIndex; ///< index of hit in the external branch
+ int IntIndex; ///< index of hit in the internal array
int iStation; ///< index of station in active stations array
double x; ///< x coordinate of position [cm]
double y; ///< y coordinate of position [cm]
diff --git a/reco/L1/CbmL1MCTrack.cxx b/reco/L1/CbmL1MCTrack.cxx
index e3330b8a10..cdf746c11e 100644
--- a/reco/L1/CbmL1MCTrack.cxx
+++ b/reco/L1/CbmL1MCTrack.cxx
@@ -23,6 +23,9 @@
#include "CbmL1.h"
#include "CbmL1Constants.h"
+#include <iomanip>
+#include <sstream>
+
#include "L1Algo/L1Algo.h"
#include "L1Algo/L1Hit.h"
@@ -193,3 +196,67 @@ void CbmL1MCTrack::CalculateIsReconstructable()
else
isAdditional = 0;
}; // bool CbmL1MCTrack::IsReconstructable()
+
+std::string CbmL1MCTrack::ToString(int verbose, bool header) const
+{
+ using std::setfill;
+ using std::setw;
+ std::stringstream msg;
+ if (header) {
+ if (verbose > 0) {
+ msg << setw(8) << setfill(' ') << "ID" << ' ';
+ msg << setw(8) << setfill(' ') << "Mother" << ' ';
+ msg << setw(8) << setfill(' ') << "PDG" << ' ';
+ if (verbose > 1) {
+ msg << setw(8) << setfill(' ') << "N h." << ' ';
+ msg << setw(8) << setfill(' ') << "N p." << ' ';
+ msg << setw(8) << setfill(' ') << "N r.tr." << ' ';
+ if (verbose > 2) { msg << setw(8) << setfill(' ') << "N t.tr." << ' '; }
+ msg << setw(12) << setfill(' ') << "zVTX [cm]" << ' ';
+ msg << setw(12) << setfill(' ') << "t [ns]" << ' ';
+ msg << setw(12) << setfill(' ') << "p [GeV/c]" << ' ';
+ }
+ msg << setw(8) << setfill(' ') << "rec-able" << ' ';
+ msg << setw(8) << setfill(' ') << "rec-ed" << ' ';
+ }
+ }
+ else {
+ if (verbose > 0) {
+ msg << setw(8) << setfill(' ') << ID << ' ';
+ msg << setw(8) << setfill(' ') << mother_ID << ' ';
+ msg << setw(8) << setfill(' ') << pdg << ' ';
+ if (verbose > 1) {
+ msg << setw(8) << setfill(' ') << Hits.size() << ' ';
+ msg << setw(8) << setfill(' ') << Points.size() << ' ';
+ msg << setw(8) << setfill(' ') << rTracks.size() << ' ';
+ if (verbose > 2) { msg << setw(8) << setfill(' ') << tTracks.size() << ' '; }
+ msg << setw(12) << setfill(' ') << z << ' ';
+ msg << setw(12) << setfill(' ') << time << ' ';
+ msg << setw(12) << setfill(' ') << p << ' ';
+ }
+ msg << setw(8) << setfill(' ') << IsReconstructable() << ' ';
+ msg << setw(8) << setfill(' ') << IsReconstructed() << ' ';
+ if (verbose > 1) {
+ msg << "\n\t- point indexes: ";
+ for (int index : Points) {
+ msg << index << ' ';
+ }
+ msg << "\n\t- hit indexes: ";
+ for (int index : Hits) {
+ msg << index << ' ';
+ }
+ msg << "\n\t- reconstructed track indexes: ";
+ for (auto* index : rTracks) {
+ msg << index << ' ';
+ }
+ if (verbose > 2) {
+ msg << "\n\t- touch track indexes: ";
+ for (auto* index : tTracks) {
+ msg << index << ' ';
+ }
+ }
+ }
+ }
+ }
+ return msg.str();
+}
diff --git a/reco/L1/CbmL1MCTrack.h b/reco/L1/CbmL1MCTrack.h
index 530fa73d2f..b7393a6eed 100644
--- a/reco/L1/CbmL1MCTrack.h
+++ b/reco/L1/CbmL1MCTrack.h
@@ -27,6 +27,7 @@
#include "TVector3.h"
#include <iostream>
+#include <string>
#include "L1Vector.h"
@@ -70,6 +71,9 @@ public:
double pt() { return sqrt(px * px + py * py); }
+ /// @brief String representation of the contents
+ std::string ToString(int verbose = 10, bool header = false) const;
+
private:
void CalculateMCCont();
void CountHitStations();
diff --git a/reco/L1/CbmL1Track.cxx b/reco/L1/CbmL1Track.cxx
new file mode 100644
index 0000000000..271fd8649b
--- /dev/null
+++ b/reco/L1/CbmL1Track.cxx
@@ -0,0 +1,60 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// @file CbmL1Track.cxx
+/// @brief Ca Tracking track representation in CBM (implementation)
+/// @since 04.04.2023
+/// @author Sergei Zharko <s.zharko@gsi.de>
+
+#include "CbmL1Track.h"
+
+#include <sstream>
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+std::string CbmL1Track::ToString(int verbose, bool header) const
+{
+ using std::setfill;
+ using std::setw;
+ std::stringstream msg;
+ if (header) {
+ msg << setw(12) << setfill(' ') << "p [GeV/c]" << ' ';
+ msg << setw(12) << setfill(' ') << "Tx" << ' ';
+ msg << setw(8) << setfill(' ') << "N hits" << ' ';
+ msg << setw(8) << setfill(' ') << "N MC tracks" << ' ';
+ msg << setw(12) << setfill(' ') << "MC tr (old)" << ' ';
+ msg << setw(12) << setfill(' ') << "MC tr (new)" << ' ';
+ msg << setw(8) << setfill(' ') << "IsGhost" << ' ';
+ msg << setw(8) << setfill(' ') << "Max.pur." << ' ';
+ msg << '\n';
+ }
+ else {
+
+
+ msg << setw(12) << setfill(' ') << GetP() << ' ';
+ msg << setw(12) << setfill(' ') << GetTx() << ' ';
+ msg << setw(8) << setfill(' ') << GetNofHits() << ' ';
+ msg << setw(8) << setfill(' ') << GetNofMCTracks() << ' ';
+ msg << setw(12) << setfill(' ') << GetMatchedMCTrackIndex() << ' ';
+ msg << setw(12) << setfill(' ') << (mcTracks.size() ? mcTracks[0]->ID : -1) << ' ';
+ msg << setw(8) << setfill(' ') << IsGhost() << ' ';
+ msg << setw(8) << setfill(' ') << GetMaxPurity() << ' ';
+ msg << '\n';
+ if (verbose > 0) {
+ msg << "\thits: ";
+ for (int iH : GetHitIndexes()) {
+ msg << iH << ' ';
+ }
+ msg << "\n\tmc tracks: (NEW)";
+ for (int iT : GetMCTrackIndexes()) {
+ msg << iT << ' ';
+ }
+ msg << "\n\tmc tracks: (Old)";
+ for (auto* iT : mcTracks) {
+ msg << iT->ID << ' ';
+ }
+ }
+ }
+ return msg.str();
+}
diff --git a/reco/L1/CbmL1Track.h b/reco/L1/CbmL1Track.h
index 6ea035c7ed..d89fd5cedc 100644
--- a/reco/L1/CbmL1Track.h
+++ b/reco/L1/CbmL1Track.h
@@ -28,7 +28,9 @@
#include "TMath.h"
+#include <iterator>
#include <map>
+#include <string>
#include <vector>
using std::map;
using std::vector;
@@ -52,13 +54,24 @@ public:
fvMcTrackIndexes.clear();
}
+ int GetNMCTracks() { return mcTracks.size(); }
+
+ bool IsGhost() const { return (maxPurity < CbmL1Constants::MinPurity); }
+
+ /// Gets charge
+ int GetCharge() const { return (T[4] > 0) ? 1 : -1; }
+
/// Gets Chi-square of track fit model
double GetChiSq() const { return chi2; }
- /// Gets a reference to hit indexes array
- const vector<int>& GetHitIndexes() const { return Hits; }
+ /// @brief Gets first hit index
+ int GetFirstHitIndex() const { return Hits.front(); }
- bool IsGhost() const { return (maxPurity < CbmL1Constants::MinPurity); }
+ /// @brief Gets last hit index
+ int GetLastHitIndex() const { return Hits.back(); }
+
+ /// Gets hit indexes
+ const auto& GetHitIndexes() const { return Hits; }
/// Gets NDF of track fit model
int GetNDF() const { return NDF; }
@@ -80,6 +93,11 @@ public:
/// Gets pointer of the associated MC track
CbmL1MCTrack* GetMCTrack() { return mcTracks[0]; }
+ /// @brief Gets index of matched MC track
+ /// @note If two tracks are matched (should not happen, if purity > 50%), the first
+ ///
+ int GetMatchedMCTrackIndex() const { return fvMcTrackIndexes.size() ? fvMcTrackIndexes[0] : -1; }
+
/// Gets number of hits of the track
int GetNofHits() const { return Hits.size(); }
@@ -90,6 +108,7 @@ public:
int GetNofStations() const { return nStations; }
/// Gets absolute value of momentum divided by the charge (absolute value) of particle [GeV/ec]
+ // NOTE: 1.e-10 precision is used in the old performance, but in FairTrackParam the 1.e-4 is used
double GetP() const { return fabs(T[4]) > 1.e-10 ? 1. / fabs(T[4]) : 1.e-10; }
/// Gets transverse momentum
@@ -110,6 +129,9 @@ public:
/// Gets z-component of momentum divided by the charge (absolute value) of particle [GeV/ec]
double GetPz() const { return std::sqrt(GetP() * GetP() / (GetTx() * GetTx() + GetTy() * GetTy() + 1)); }
+ /// Gets track polar angle
+ double GetTheta() const { return std::acos(1. / std::sqrt(GetTx() * GetTx() + GetTy() * GetTy() + 1)); }
+
/// Gets track slope along x-axis
double GetTx() const { return T[2]; }
@@ -126,6 +148,10 @@ public:
static bool comparePChi2(const CbmL1Track* a, const CbmL1Track* b) { return (a->chi2 < b->chi2); }
+ /// @brief Provides a string representation of object
+ /// @param verbose Verbosity level
+ /// @param header If true, header will be printed
+ std::string ToString(int verbose = 10, bool header = false) const;
double Tpv[L1TrackPar::kNparTr]; ///< Track parameters at primary vertex
double Cpv[L1TrackPar::kNparCov]; ///< Track covariance matrix at primary vertex
diff --git a/reco/L1/L1Algo/L1Constants.h b/reco/L1/L1Algo/L1Constants.h
index af74194f3a..08047172e5 100644
--- a/reco/L1/L1Algo/L1Constants.h
+++ b/reco/L1/L1Algo/L1Constants.h
@@ -41,6 +41,8 @@ namespace L1Constants
constexpr int kMaxNthreads = 1u << kThreadBits; ///< Max number of threads, 2^6 = 64
constexpr int kMaxNtriplets = 1u << kTripletBits; ///< Max number of triplets, 2^20 = 1,048,576
+ constexpr uint8_t kDetBits = 4u; // Max 16 detectors
+
/// Max number of track groups
/// NOTE: For a "track group" definition see L1Parameters.h, GetSearchWindow function
constexpr int kMaxNtrackGroups = 4;
@@ -70,28 +72,36 @@ namespace L1Constants
} // namespace control
- /// Physics constants
+ // Physics constants
namespace phys
{
/* Particle masses used for track fit */
- constexpr float kMuonMass = 0.10565800f; ///< Muon mass [GeV/c2]
- constexpr float kElectronMass = 0.000511f; ///< Electron mass [GeV/c2]
+ constexpr float kMuonMass = 0.10565800f; ///< Muon mass [GeV/c2] !! TODO: discrepancy in last two digits
+ constexpr float kPionMass = 0.13957039f; ///< Pion mass [GeV/c2]
+ constexpr float kElectronMass = 0.000511f; ///< Electron mass [GeV/c2]
+ constexpr float kProtonMass = 0.93827209f; ///< Proton mass [GeV/c2] (0.93827208816 - PDG 11.08.2022)
constexpr double kSpeedOfLight = 29.9792458; ///< Speed of light [cm/ns]
- } // namespace phys
+ constexpr float kSpeedOfLightF = 29.9792; ///< Speed of light [cm/ns] (single precision)
+ constexpr double kSpeedOfLightInv = 1 / kSpeedOfLight; ///< Inverse speed of light [cm/ns]
+ constexpr float kSpeedOfLightInvF = 1 / kSpeedOfLightF; ///< Inverse speed of light [cm/ns] (single precision)
+ } // namespace phys
- /// Math constants
+ // Math constants
namespace math
{
constexpr double kPi = 3.14159265358979323846; ///< Value of PI, used in ROOT TMath
}
- /// Miscellaneous constants
+ // Miscellaneous constants
namespace misc
{
constexpr int kAssertionLevel = 0; ///< Assertion level
constexpr int kAlignment = 16; ///< Default alignment of data (bytes)
} // namespace misc
+ // Units
+
+
// Colors of terminal log
namespace clrs
{
diff --git a/reco/L1/L1Algo/L1Field.cxx b/reco/L1/L1Algo/L1Field.cxx
index c91cf325b3..2dbaba5840 100644
--- a/reco/L1/L1Algo/L1Field.cxx
+++ b/reco/L1/L1Algo/L1Field.cxx
@@ -182,7 +182,9 @@ void L1FieldRegion::Get(const fvec& x, const fvec& y, const fvec& z, fvec* B) co
for (size_t i = 0; i < fvec::size(); i++) {
double inPos[3] = {x[i], y[i], z[i]};
double outB[3];
- CbmKF::Instance()->GetMagneticField()->GetFieldValue(inPos, outB);
+ assert(CbmKF::Instance());
+ assert(CbmKF::Instance()->GetMagneticField());
+ CbmKF::Instance()->GetMagneticField()->GetFieldValue(inPos, outB); // TODO: replace with functional object
B[0][i] = outB[0];
B[1][i] = outB[1];
B[2][i] = outB[2];
diff --git a/reco/L1/L1Algo/L1Fit.h b/reco/L1/L1Algo/L1Fit.h
index 458d9d576a..2be29385f5 100644
--- a/reco/L1/L1Algo/L1Fit.h
+++ b/reco/L1/L1Algo/L1Fit.h
@@ -41,7 +41,7 @@ public:
}
template<typename T>
- void SetTrack(const T p[7], const T c[28])
+ void SetTrack(const T p[8], const T c[28])
{
fTr.copyFromArrays(p, c);
fQp0 = fTr.qp;
diff --git a/reco/L1/L1Algo/L1TrackPar.h b/reco/L1/L1Algo/L1TrackPar.h
index 6fe8491644..a8a118b67a 100644
--- a/reco/L1/L1Algo/L1TrackPar.h
+++ b/reco/L1/L1Algo/L1TrackPar.h
@@ -1,12 +1,13 @@
-/* Copyright (C) 2007-2020 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+/* Copyright (C) 2007-2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
SPDX-License-Identifier: GPL-3.0-only
- Authors: Igor Kulakov [committer], Maksym Zyzak */
+ Authors: Igor Kulakov [committer], Maksym Zyzak, Sergei Zharko */
#ifndef L1TrackPar_h
#define L1TrackPar_h 1
#include "CaSimd.h"
#include "L1Def.h"
+#include "L1Undef.h"
using namespace cbm::algo::ca;
@@ -41,7 +42,8 @@ public:
template<typename T>
L1TrackPar(const T* tr, const T* C)
{
- Set(tr, C);
+ //Set(tr, C); // Not implemented
+ copyFromArrays(tr, C);
}
//template<typename T>
@@ -63,8 +65,62 @@ public:
return c[ind];
}
+ // ** Parameter getters **
+
+ /// @brief Gets x position [cm]
+ fscal GetX() const { return x[0]; }
+
+ /// @brief Gets x position error [cm]
+ fscal GetXErr() const { return (std::isfinite(C00[0]) && C00[0] > 0) ? std::sqrt(C00[0]) : undef::kF; }
+
+ /// @brief Gets y position [cm]
+ fscal GetY() const { return y[0]; }
+
+ /// @brief Gets y position error [cm]
+ fscal GetYErr() const { return (std::isfinite(C11[0]) && C11[0] > 0) ? std::sqrt(C11[0]) : undef::kF; }
+
+ /// @brief Gets slope along x-axis
+ fscal GetTx() const { return tx[0]; }
+
+ /// @brief Gets error of slope along x-axis
+ fscal GetTxErr() const { return (std::isfinite(C22[0]) && C22[0] > 0) ? std::sqrt(C22[0]) : undef::kF; }
+
+ /// @brief Gets slope along y-axis
+ fscal GetTy() const { return ty[0]; }
+
+ /// @brief Gets error of slope along y-axis
+ fscal GetTyErr() const { return (std::isfinite(C33[0]) && C33[0] > 0) ? std::sqrt(C33[0]) : undef::kF; }
+
+ /// @brief Gets charge over momentum [ec/GeV]
+ fscal GetQp() const { return qp[0]; }
+
+ /// @brief Gets error of charge over momentum [ec/GeV]
+ fscal GetQpErr() const { return (std::isfinite(C44[0]) && C44[0] > 0) ? std::sqrt(C44[0]) : undef::kF; }
+
+ /// @brief Gets time [ns]
+ fscal GetTime() const { return t[0]; }
+
+ /// @brief Gets time error [ns]
+ fscal GetTimeErr() const { return (std::isfinite(C55[0]) && C55[0] > 0) ? std::sqrt(C55[0]) : undef::kF; }
+
+ /// @brief Gets inverse speed [ns/cm]
+ fscal GetInvSpeed() const { return vi[0]; }
+
+ /// @brief Gets inverse speed error [ns/cm]
+ fscal GetInvSpeedErr() const { return (std::isfinite(C66[0]) && C66[0] > 0) ? std::sqrt(C66[0]) : undef::kF; }
+
+ /// @brief Resets variances of track parameters
+ /// @param c00 Variance of x-position [cm2]
+ /// @param c11 Variance of y-position [cm2]
+ /// @param c22 Variance of slope along x-axis
+ /// @param c33 Variance of slope along y-axis
+ /// @param c44 Variance of charge over momentum [(ec/GeV)2]
+ /// @param c55 Variance of time [ns2]
+ /// @param c66 Variance of inverse speed [1/c2]
void ResetErrors(fvec c00, fvec c11, fvec c22, fvec c33, fvec c44, fvec c55, fvec c66);
+ /// @brief Prints parameters
+ /// @param i Index of SIMD vector element (if -1, the whole vector is printed out)
void Print(int i = -1) const;
void PrintCorrelations(int i = -1) const;
diff --git a/reco/L1/catools/CaToolsMCPoint.h b/reco/L1/catools/CaToolsMCPoint.h
index 83aec7054b..d1476138d4 100644
--- a/reco/L1/catools/CaToolsMCPoint.h
+++ b/reco/L1/catools/CaToolsMCPoint.h
@@ -13,38 +13,39 @@
#include <string>
#include "CaToolsLinkKey.h"
+#include "L1Constants.h"
#include "L1Undef.h"
#include "L1Vector.h"
-
enum class L1DetectorID;
-/// Class describes a Monte-Carlo point used in CA tracking QA analysis
namespace ca::tools
{
+ /// @brief Class describes a unified MC-point, used in CA tracking QA analysis
+ ///
class MCPoint {
public:
- /// Default constructor
+ /// @brief Default constructor
MCPoint() = default;
- /// Destructor
+ /// @brief Destructor
~MCPoint() = default;
- /// Copy constructor
+ /// @brief Copy constructor
MCPoint(const MCPoint&) = default;
- /// Move constructor
+ /// @brief Move constructor
MCPoint(MCPoint&&) = default;
- /// Copy assignment operator
+ /// @brief Copy assignment operator
MCPoint& operator=(const MCPoint&) = default;
- /// Move assignment operator
+ /// @brief Move assignment operator
MCPoint& operator=(MCPoint&&) = default;
- /// Adds index of hits from the container of hits of event/TS
- /// \param iH A hit index in the external hits container of event/TS
+ /// @brief Adds index of hits from the container of hits of event/TS
+ /// @param iH A hit index in the external hits container of event/TS
void AddHitID(int iH) { fvHitIndexes.push_back_no_warning(iH); }
@@ -52,211 +53,257 @@ namespace ca::tools
// ** Getters **
// *********************
- /// Gets charge of the particle [e]
+ /// @brief Gets charge of the particle [e]
double GetCharge() const { return fCharge; }
- /// Gets detector ID
+ /// @brief Gets detector ID
L1DetectorID GetDetectorId() const { return fDetectorId; }
- /// Gets MC event ID
+ /// @brief Gets MC event ID
int GetEventId() const { return fLinkKey.fEvent; }
- /// Gets MC file ID
+ /// @brief Gets MC file ID
int GetFileId() const { return fLinkKey.fFile; }
- /// Gets index of this point in internal CA container
+ /// @brief Gets index of this point in internal CA container
int GetId() const { return fId; }
- /// Gets container of matched hit indexes
+ /// @brief Gets inverse speed at reference z of station [ns/cm]
+ int GetInvSpeed() const
+ {
+ return std::sqrt(1. + GetMass() * GetMass() / GetP() / GetP()) * L1Constants::phys::kSpeedOfLightInv;
+ }
+
+ /// @brief Gets inverse speed at entrance to station [ns/cm]
+ int GetInvSpeedIn() const
+ {
+ return std::sqrt(1. + GetMass() * GetMass() / GetPIn() / GetPIn()) * L1Constants::phys::kSpeedOfLightInv;
+ }
+
+ /// @brief Gets inverse speed at exit of station [ns/cm]
+ int GetInvSpeedOut() const
+ {
+ return std::sqrt(1. + GetMass() * GetMass() / GetPOut() / GetPOut()) * L1Constants::phys::kSpeedOfLightInv;
+ }
+
+ /// @brief Gets container of matched hit indexes
const auto& GetHitIndexes() const { return fvHitIndexes; }
- /// Gets link key
+ /// @brief Gets link key
LinkKey GetLinkKey() const { return fLinkKey; }
- /// Gets mass of the particle [GeV/c2]
+ /// @brief Gets mass of the particle [GeV/c2]
double GetMass() const { return fMass; }
- /// Gets mother ID of the track
+ /// @brief Gets mother ID of the track
int GetMotherId() const { return fMotherId; }
- /// Gets track momentum absolute value at reference z of station [GeV/c]
+ /// @brief Gets track momentum absolute value at reference z of station [GeV/c]
double GetP() const { return std::sqrt(fMom[0] * fMom[0] + fMom[1] * fMom[1] + fMom[2] * fMom[2]); }
- /// Gets PDG code of the particle
+ /// @brief Gets PDG code of the particle
int GetPdgCode() const { return fPdgCode; }
- /// Gets track momentum absolute value at entrance to station [GeV/c]
+ /// @brief Gets track momentum absolute value at entrance to station [GeV/c]
double GetPIn() const { return std::sqrt(fMomIn[0] * fMomIn[0] + fMomIn[1] * fMomIn[1] + fMomIn[2] * fMomIn[2]); }
- /// Gets track momentum absolute value at exit of station [GeV/c]
+ /// @brief Gets track momentum absolute value at exit of station [GeV/c]
double GetPOut() const
{
return std::sqrt(fMomOut[0] * fMomOut[0] + fMomOut[1] * fMomOut[1] + fMomOut[2] * fMomOut[2]);
}
- /// Gets track momentum x component at reference z of station [GeV/c]
+ /// @brief Gets track momentum x component at reference z of station [GeV/c]
double GetPx() const { return fMom[0]; }
- /// Gets track momentum x component at entrance to station [GeV/c]
+ /// @brief Gets track momentum x component at entrance to station [GeV/c]
double GetPxIn() const { return fMomIn[0]; }
- /// Gets track momentum x component at exit of station [GeV/c]
+ /// @brief Gets track momentum x component at exit of station [GeV/c]
double GetPxOut() const { return fMomOut[0]; }
- /// Gets track momentum y component at reference z of station [GeV/c]
+ /// @brief Gets track momentum y component at reference z of station [GeV/c]
double GetPy() const { return fMom[1]; }
- /// Gets track momentum y component at entrance to station [GeV/c]
+ /// @brief Gets track momentum y component at entrance to station [GeV/c]
double GetPyIn() const { return fMomIn[1]; }
- /// Gets track momentum y component at exit of station [GeV/c]
+ /// @brief Gets track momentum y component at exit of station [GeV/c]
double GetPyOut() const { return fMomOut[1]; }
- /// Gets track momentum z component at reference z of station [GeV/c]
+ /// @brief Gets track momentum z component at reference z of station [GeV/c]
double GetPz() const { return fMom[2]; }
- /// Gets track momentum z component at entrance to station [GeV/c]
+ /// @brief Gets track momentum z component at entrance to station [GeV/c]
double GetPzIn() const { return fMomIn[2]; }
- /// Gets track momentum z component at exit of station [GeV/c]
+ /// @brief Gets track momentum z component at exit of station [GeV/c]
double GetPzOut() const { return fMomOut[2]; }
- /// Gets global ID of the active tracking station
+ /// @brief Gets track charge over momentum at reference z of station [ec/GeV]
+ double GetQp() const { return fCharge / GetP(); }
+
+ /// @brief Gets track momentum absolute value at entrance to station [ec/GeV]
+ double GetQpIn() const { return fCharge / GetPIn(); }
+
+ /// @brief Gets track momentum absolute value at exit of station [ec/GeV]
+ double GetQpOut() const { return fCharge / GetPOut(); }
+
+ /// @brief Gets global ID of the active tracking station
int GetStationId() const { return fStationId; }
- /// Gets time [ns]
+ /// @brief Gets time [ns]
double GetTime() const { return fTime; }
- /// Gets ID of track from the internal CA MC track container (within event/TS)
+ /// @brief Gets ID of track from the internal CA MC track container (within event/TS)
int GetTrackId() const { return fTrackId; }
- /// Gets x coordinate at reference z of station [cm]
+ /// @brief Gets slope along x-axis at reference z of station
+ double GetTx() const { return fMom[0] / fMom[2]; }
+
+ /// @brief Gets slope along x-axis at entrance to station
+ double GetTxIn() const { return fMomIn[0] / fMomIn[2]; }
+
+ /// @brief Gets slope along x-axis at exit of station
+ double GetTxOut() const { return fMomOut[0] / fMomOut[2]; }
+
+ /// @brief Gets slope along x-axis at reference z of station
+ double GetTy() const { return fMom[1] / fMom[2]; }
+
+ /// @brief Gets slope along x-axis at entrance to station
+ double GetTyIn() const { return fMomIn[1] / fMomIn[2]; }
+
+ /// @brief Gets slope along x-axis at exit of station
+ double GetTyOut() const { return fMomOut[1] / fMomOut[2]; }
+
+ /// @brief Gets x coordinate at reference z of station [cm]
double GetX() const { return fPos[0]; }
- /// Gets x coordinate at entrance to station [cm]
+ /// @brief Gets x coordinate at entrance to station [cm]
double GetXIn() const { return fPosIn[0]; }
- /// Gets x coordinate at exit of station [cm]
+ /// @brief Gets x coordinate at exit of station [cm]
double GetXOut() const { return fPosOut[0]; }
- /// Gets y coordinate at reference z of station [cm]
+ /// @brief Gets y coordinate at reference z of station [cm]
double GetY() const { return fPos[1]; }
- /// Gets y coordinate at entrance to station [cm]
+ /// @brief Gets y coordinate at entrance to station [cm]
double GetYIn() const { return fPosIn[1]; }
- /// Gets y coordinate at exit of station [cm]
+ /// @brief Gets y coordinate at exit of station [cm]
double GetYOut() const { return fPosOut[1]; }
- /// Gets z coordinate at reference z of station [cm]
+ /// @brief Gets z coordinate at reference z of station [cm]
double GetZ() const { return fPos[2]; }
- /// Gets z coordinate at entrance to station [cm]
+ /// @brief Gets z coordinate at entrance to station [cm]
double GetZIn() const { return fPosIn[2]; }
- /// Gets z coordinate at exit of station [cm]
+ /// @brief Gets z coordinate at exit of station [cm]
double GetZOut() const { return fPosOut[2]; }
// *********************
// ** Setters **
// *********************
- /// Sets particle charge [e]
+ /// @brief Sets particle charge [e]
void SetCharge(double charge) { fCharge = charge; }
- /// Sets detector ID
+ /// @brief Sets detector ID
void SetDetectorId(L1DetectorID detId) { fDetectorId = detId; }
- /// Sets index of MC event containing this point
+ /// @brief Sets index of MC event containing this point
void SetEventId(int eventId) { fLinkKey.fEvent = eventId; }
- /// Sets index of this point in external data structures
+ /// @brief Sets index of this point in external data structures
void SetExternalId(int id) { fLinkKey.fIndex = id; }
- /// Sets index of MC file containing this point
+ /// @brief Sets index of MC file containing this point
void SetFileId(int fileId) { fLinkKey.fFile = fileId; }
- /// Sets index of this point in the CA internal structure
+ /// @brief Sets index of this point in the CA internal structure
void SetId(int id) { fId = id; }
- /// Sets particle mass [GeV/c2]
+ /// @brief Sets particle mass [GeV/c2]
void SetMass(double mass) { fMass = mass; }
- /// Sets index of mother track in the internal CA data structures
+ /// @brief Sets index of mother track in the internal CA data structures
void SetMotherId(int motherId) { fMotherId = motherId; }
- /// Sets PDG code
+ /// @brief Sets PDG code
void SetPdgCode(int pdg) { fPdgCode = pdg; }
- /// Sets track momentum x component at reference z of station [GeV/c]
+ /// @brief Sets track momentum x component at reference z of station [GeV/c]
void SetPx(double px) { fMom[0] = px; }
- /// Sets track momentum x component at entrance to station [GeV/c]
+ /// @brief Sets track momentum x component at entrance to station [GeV/c]
void SetPxIn(double px) { fMomIn[0] = px; }
- /// Sets track momentum x component at exit of station [GeV/c]
+ /// @brief Sets track momentum x component at exit of station [GeV/c]
void SetPxOut(double px) { fMomOut[0] = px; }
- /// Sets track momentum y component at reference z of station [GeV/c]
+ /// @brief Sets track momentum y component at reference z of station [GeV/c]
void SetPy(double py) { fMom[1] = py; }
- /// Sets track momentum y component at entrance to station [GeV/c]
+ /// @brief Sets track momentum y component at entrance to station [GeV/c]
void SetPyIn(double py) { fMomIn[1] = py; }
- /// Sets track momentum y component at exit of station [GeV/c]
+ /// @brief Sets track momentum y component at exit of station [GeV/c]
void SetPyOut(double py) { fMomOut[1] = py; }
- /// Sets track momentum z component at reference z of station [GeV/c]
+ /// @brief Sets track momentum z component at reference z of station [GeV/c]
void SetPz(double pz) { fMom[2] = pz; }
- /// Sets track momentum z component at entrance to station [GeV/c]
+ /// @brief Sets track momentum z component at entrance to station [GeV/c]
void SetPzIn(double pz) { fMomIn[2] = pz; }
- /// Sets track momentum z component at exit of station [GeV/c]
+ /// @brief Sets track momentum z component at exit of station [GeV/c]
void SetPzOut(double pz) { fMomOut[2] = pz; }
- /// Sets global index of active station
+
+ /// @brief Sets global index of active station
void SetStationId(int stationId) { fStationId = stationId; }
- /// Sets time [ns]
+ /// @brief Sets time [ns]
void SetTime(double time) { fTime = time; }
- /// Sets track ID in the CA internal track container (within event/TS)
+ /// @brief Sets track ID in the CA internal track container (within event/TS)
void SetTrackId(int trackId) { fTrackId = trackId; }
- /// Sets x coordinate at reference z of station [cm]
+ /// @brief Sets x coordinate at reference z of station [cm]
void SetX(double x) { fPos[0] = x; }
- /// Sets x coordinate at entrance to station [cm]
+ /// @brief Sets x coordinate at entrance to station [cm]
void SetXIn(double x) { fPosIn[0] = x; }
- /// Sets x coordinate at exit of station [cm]
+ /// @brief Sets x coordinate at exit of station [cm]
void SetXOut(double x) { fPosOut[0] = x; }
- /// Sets y coordinate at reference z of station [cm]
+ /// @brief Sets y coordinate at reference z of station [cm]
void SetY(double y) { fPos[1] = y; }
- /// Sets y coordinate at entrance to station [cm]
+ /// @brief Sets y coordinate at entrance to station [cm]
void SetYIn(double y) { fPosIn[1] = y; }
- /// Sets x coordinate at exit of station [cm]
+ /// @brief Sets x coordinate at exit of station [cm]
void SetYOut(double y) { fPosOut[1] = y; }
- /// Sets z coordinate at reference z of station [cm]
+ /// @brief Sets z coordinate at reference z of station [cm]
void SetZ(double z) { fPos[2] = z; }
- /// Sets z coordinate at entrance to station [cm]
+ /// @brief Sets z coordinate at entrance to station [cm]
void SetZIn(double z) { fPosIn[2] = z; }
- /// Sets z coordinate at exit of station [cm]
+ /// @brief Sets z coordinate at exit of station [cm]
void SetZOut(double z) { fPosOut[2] = z; }
- /// Prints content for a given verbosity level
- /// \param verbose Verbosity level:
+ /// @brief Prints content for a given verbosity level
+ /// @param verbose Verbosity level:
/// -#0: Prints nothing
/// -#1: Prints track ID, station ID, time and position
/// -#2: Also prints zIn, zOut, absolute momentum, as well as point, event and file IDs
- /// \param printHeader If true, parameter names will be printed instead of the parameters themselves
+ /// @param printHeader If true, parameter names will be printed instead of the parameters themselves
std::string ToString(int verbose, bool printHeader = false) const;
// ****************************
@@ -285,7 +332,7 @@ namespace ca::tools
L1DetectorID fDetectorId; ///< Detector ID of MC point
- L1Vector<int> fvHitIndexes {"ca::tools::MCPoint::fvHitIndexes"};
+ L1Vector<int> fvHitIndexes {"ca::tools::MCPoint::fvHitIndexes"}; ///< Indexes of hits, assigned to this point
};
} // namespace ca::tools
diff --git a/reco/L1/qa/CbmCaOutputQa.cxx b/reco/L1/qa/CbmCaOutputQa.cxx
index 51080c550a..b0c7ba4004 100644
--- a/reco/L1/qa/CbmCaOutputQa.cxx
+++ b/reco/L1/qa/CbmCaOutputQa.cxx
@@ -36,7 +36,8 @@ void OutputQa::FillHistograms()
// ** Fill distributions for reconstructed tracks (including ghost ones) **
// ************************************************************************
- for (const auto& recoTrk : fvRecoTracks) {
+ for (size_t iTrkReco = 0; iTrkReco < fvRecoTracks.size(); ++iTrkReco) {
+ const auto& recoTrk = fvRecoTracks[iTrkReco];
// Reject tracks, which do not contain hits
if (recoTrk.GetNofHits() < 1) { continue; }
@@ -84,40 +85,90 @@ void OutputQa::FillHistograms()
}
} // recoTrk.GetNDF() > 0: end
- if (recoTrk.IsGhost()) {
- fph_ghost_purity->Fill(100 * recoTrk.GetMaxPurity());
- fph_ghost_p->Fill(recoTrk.GetP());
- fph_ghost_pt->Fill(recoTrk.GetPt());
- fph_ghost_phi->Fill(recoTrk.GetPhi());
- fph_ghost_nhits->Fill(recoTrk.GetNofHits());
- fph_ghost_fsta->Fill(hitFst.GetStationId());
- fph_ghost_fhitR->Fill(hitFst.GetR());
-
- // z-positions of the first and second hit stations
- double z0 = fpParameters->GetStation(hitFst.GetStationId()).fZ[0];
- double z1 = fpParameters->GetStation(hitSnd.GetStationId()).fZ[0];
-
- if (fabs(z1 - z0) > 1.e-4) { // test for z1 != z2
- fph_ghost_tx->Fill((hitSnd.GetX() - hitFst.GetX()) / (z1 - z0));
- fph_ghost_ty->Fill((hitSnd.GetY() - hitFst.GetY()) / (z1 - z0));
- }
+ if (IsMCUsed()) {
+ if (recoTrk.IsGhost()) {
+ fph_ghost_purity->Fill(100 * recoTrk.GetMaxPurity());
+ fph_ghost_p->Fill(recoTrk.GetP());
+ fph_ghost_pt->Fill(recoTrk.GetPt());
+ fph_ghost_phi->Fill(recoTrk.GetPhi());
+ fph_ghost_nhits->Fill(recoTrk.GetNofHits());
+ fph_ghost_fsta->Fill(hitFst.GetStationId());
+ fph_ghost_fhitR->Fill(hitFst.GetR());
+
+ // z-positions of the first and second hit stations
+ double z0 = fpParameters->GetStation(hitFst.GetStationId()).fZ[0];
+ double z1 = fpParameters->GetStation(hitSnd.GetStationId()).fZ[0];
+
+ if (fabs(z1 - z0) > 1.e-4) { // test for z1 != z2
+ fph_ghost_tx->Fill((hitSnd.GetX() - hitFst.GetX()) / (z1 - z0));
+ fph_ghost_ty->Fill((hitSnd.GetY() - hitFst.GetY()) / (z1 - z0));
+ }
+
+ fph_ghost_nhits_vs_p->Fill(recoTrk.GetP(), recoTrk.GetNofHits());
+ fph_ghost_fsta_vs_p->Fill(recoTrk.GetP(), hitFst.GetStationId());
+ if (hitFst.GetStationId() <= hitLst.GetStationId()) {
+ fph_ghost_lsta_vs_fsta->Fill(hitFst.GetStationId(), hitLst.GetStationId());
+ }
+ } // recoTrk.IsGhost(): end
+ }
+
+ // ********************************
+ // ** Different track categories **
+ // ********************************
+
+ fvpTrackHistograms[ETrackType::kAll]->FillRecoTrack(iTrkReco);
+
+ if (IsMCUsed()) {
+ int iTrkMC = recoTrk.GetMatchedMCTrackIndex();
+ if (iTrkMC > -1) {
+ const auto& mcTrk = fMCData.GetTrack(iTrkMC);
+ int pdg = mcTrk.GetPdgCode();
+ bool isPrimary = mcTrk.IsPrimary();
- fph_ghost_nhits_vs_p->Fill(recoTrk.GetP(), recoTrk.GetNofHits());
- fph_ghost_fsta_vs_p->Fill(recoTrk.GetP(), hitFst.GetStationId());
- if (hitFst.GetStationId() <= hitLst.GetStationId()) {
- fph_ghost_lsta_vs_fsta->Fill(hitFst.GetStationId(), hitLst.GetStationId());
+ if (isPrimary) {
+ fvpTrackHistograms[ETrackType::kPrim]->FillRecoTrack(iTrkReco);
+
+ if (pdg == +211) { fvpTrackHistograms[ETrackType::kPrimPIP]->FillRecoTrack(iTrkReco); }
+ if (pdg == -211) { fvpTrackHistograms[ETrackType::kPrimPIM]->FillRecoTrack(iTrkReco); }
+ }
+ else {
+ fvpTrackHistograms[ETrackType::kSec]->FillRecoTrack(iTrkReco);
+
+ if (pdg == +211) { fvpTrackHistograms[ETrackType::kSecPIP]->FillRecoTrack(iTrkReco); }
+ if (pdg == -211) { fvpTrackHistograms[ETrackType::kSecPIM]->FillRecoTrack(iTrkReco); }
+ }
}
- } // recoTrk.IsGhost(): end
- } // loop over recoTrk: end
+ }
+
+
+ } // loop over recoTrk: end
// *************************************
// ** Fill distributions of MC-tracks **
// *************************************
+ if (IsMCUsed()) {
+ for (int iTrkMC = 0; iTrkMC < fMCData.GetNofTracks(); ++iTrkMC) {
+ const auto& mcTrk = fMCData.GetTrack(iTrkMC);
+ fvpTrackHistograms[ETrackType::kAll]->FillMCTrack(iTrkMC);
+
+ int pdg = mcTrk.GetPdgCode();
+ bool isPrimary = mcTrk.IsPrimary();
+
+ if (isPrimary) {
+ fvpTrackHistograms[ETrackType::kPrim]->FillMCTrack(iTrkMC);
+
+ if (pdg == +211) { fvpTrackHistograms[ETrackType::kPrimPIP]->FillMCTrack(iTrkMC); }
+ if (pdg == -211) { fvpTrackHistograms[ETrackType::kPrimPIM]->FillMCTrack(iTrkMC); }
+ }
+ else {
+ fvpTrackHistograms[ETrackType::kSec]->FillMCTrack(iTrkMC);
- //for (const auto& mcTrk: fpMCData->GetTrackContainer()) {
- //
- //}
+ if (pdg == +211) { fvpTrackHistograms[ETrackType::kSecPIP]->FillMCTrack(iTrkMC); }
+ if (pdg == -211) { fvpTrackHistograms[ETrackType::kSecPIM]->FillMCTrack(iTrkMC); }
+ }
+ }
+ }
}
// ---------------------------------------------------------------------------------------------------------------------
@@ -199,6 +250,25 @@ InitStatus OutputQa::InitDataBranches()
//
InitStatus OutputQa::InitHistograms()
{
+ auto RegisterTrackQa = [&](const char* typeName, ETrackType type, bool bSuppressMC = false) {
+ bool bUseMC = IsMCUsed() && !bSuppressMC;
+ fvpTrackHistograms[type] = std::make_unique<TrackTypeQa>(typeName, fsPrefix.Data(), bUseMC, fpFolderRoot);
+ fvpTrackHistograms[type]->RegisterParameters(fpParameters);
+ fvpTrackHistograms[type]->RegisterRecoHits(fvHits);
+ fvpTrackHistograms[type]->RegisterRecoTracks(fvRecoTracks);
+ fvpTrackHistograms[type]->RegisterMCData(fMCData);
+ fvpTrackHistograms[type]->Init();
+ };
+
+ RegisterTrackQa("all", ETrackType::kAll);
+ RegisterTrackQa("prim", ETrackType::kPrim);
+ RegisterTrackQa("sec", ETrackType::kSec);
+ RegisterTrackQa("prim_pip", ETrackType::kPrimPIP);
+ RegisterTrackQa("prim_pim", ETrackType::kPrimPIM);
+ RegisterTrackQa("sec_pip", ETrackType::kSecPIP);
+ RegisterTrackQa("sec_pim", ETrackType::kSecPIM);
+
+
fph_reco_p = MakeHistoFromConfig<TH1F>("reco_p");
fph_reco_pt = MakeHistoFromConfig<TH1F>("reco_pt");
fph_reco_phi = MakeHistoFromConfig<TH1F>("reco_phi");
@@ -230,6 +300,11 @@ InitStatus OutputQa::InitHistograms()
fph_ghost_fsta_vs_p = MakeHistoFromConfig<TH2F>("ghost_fsta_vs_p");
fph_ghost_lsta_vs_fsta = MakeHistoFromConfig<TH2F>("ghost_lsta_vs_fsta");
+ // Reconstructed tracks vs. MC quantities
+ fph_reco_pMC = MakeHistoFromConfig<TH1F>("reco_pMC");
+ fph_reco_yMC = MakeHistoFromConfig<TH1F>("reco_yMC");
+ fph_reco_pMC_vs_yMC = MakeHistoFromConfig<TH2F>("reco_pMC_vs_yMC");
+
// Residuals and pools of track parameters
fph_fst_res_x = MakeHistoFromConfig<TH1F>("fst_res_x");
fph_fst_res_y = MakeHistoFromConfig<TH1F>("fst_res_y");
diff --git a/reco/L1/qa/CbmCaOutputQa.h b/reco/L1/qa/CbmCaOutputQa.h
index 91b54aca2d..930134f741 100644
--- a/reco/L1/qa/CbmCaOutputQa.h
+++ b/reco/L1/qa/CbmCaOutputQa.h
@@ -12,10 +12,14 @@
#include "CbmCaMCModule.h"
#include "CbmCaTimeSliceReader.h"
+#include "CbmCaTrackTypeQa.h"
#include "CbmL1DetectorID.h"
#include "CbmL1Hit.h"
#include "CbmQaTask.h"
+#include <array>
+#include <memory>
+
#include "CaToolsDebugger.h"
#include "L1Parameters.h"
@@ -25,6 +29,52 @@ namespace cbm::ca
///
class OutputQa : public CbmQaTask {
public:
+ /// @brief Enumeration of track category
+ enum ETrackType
+ {
+ kAll, ///< all tracks
+ kPrim, ///< primary tracks
+ kSec, ///< secondary tracks
+ // kPrimEP, ///< primary electron tracks
+ // kPrimEM, ///< primary positron tracks
+ // kSecEP, ///< secondary electron tracks
+ // kSecEM, ///< secondary positron tracks
+ // kPrimMUP, ///< primary .....
+ // kPrimMUM,
+ // kSecMUP,
+ // kSecMUM,
+ kPrimPIP,
+ kPrimPIM,
+ kSecPIP,
+ kSecPIM,
+ // kPrimKP,
+ // kPrimKM,
+ // kSecKP,
+ // kSecKM,
+ // kPrimP,
+ // kPrimPBAR,
+ // kSecP,
+ // kSecPBAR,
+ // kPrimD,
+ // kPrimDBAR,
+ // kSecD,
+ // kSecDBAR,
+ // kPrimT,
+ // kPrimTBAR,
+ // kSecT,
+ // kSecTBAR,
+ // kPrim3HE,
+ // kPrim3HEBAR,
+ // kSec3HE,
+ // kSec3HEBAR,
+ // kPrim4HE,
+ // kPrim4HEBAR,
+ // kSec4HE,
+ // kSec4HEBAR,
+ kEND
+ };
+
+
/// @brief Constructor from parameters
/// @param verbose Verbosity level
/// @param isMCUsed Flag, if MC information is available for this task
@@ -122,7 +172,15 @@ namespace cbm::ca
::ca::tools::MCData fMCData; ///< Input MC data (points and tracks)
- // ** Pointers to histogram objects **
+ // *************************
+ // ** List of histograms **
+ // *************************
+
+ /// Histograms of different track types
+ std::array<std::unique_ptr<TrackTypeQa>, ETrackType::kEND> fvpTrackHistograms;
+
+ std::array<std::tuple<Style_t, Marker_t, Color_t>, ETrackType::kEND> fvpTrackHistogramStyles;
+
// Reconstructed tracks
TH1F* fph_reco_p = nullptr; ///< Total momentum over charge of reconstructed tracks
TH1F* fph_reco_pt = nullptr; ///< Transverse momentum over charge of reconstructed tracks
@@ -132,7 +190,6 @@ namespace cbm::ca
TH1F* fph_reco_fhitR = nullptr; ///< Distance of the first hit from z-axis for reconstructed tracks
TH1F* fph_reco_nhits = nullptr; ///< Hit number of reconstructed tracks
TH1F* fph_reco_fsta = nullptr; ///< First station index of reconstructed tracks
- TH1F* fph_reco_purity = nullptr; ///< Purity of reconstructed tracks (\note purity requires MC information)
TH1F* fph_reco_chi2_ndf = nullptr; ///< Fit chi2 over NDF of reconstructed tracks
TH1F* fph_reco_prob = nullptr; ///< Fit probability of reconstructed tracks
TH1F* fph_rest_chi2_ndf = nullptr; ///< Fit chi2 over NDF of non-reconstructable tracks
@@ -154,6 +211,12 @@ namespace cbm::ca
TH2F* fph_ghost_fsta_vs_p = nullptr; ///< First station index vs. total momentum over charge for ghost tracks
TH2F* fph_ghost_lsta_vs_fsta = nullptr; ///< Last station index vs. first station index of ghost tracks
+ // Reco tracks vs MC variables
+ TH1F* fph_reco_pMC = nullptr; ///< Reconstructed track distribution vs MC momemntum
+ TH1F* fph_reco_yMC = nullptr; ///< Reconstructed track distribution vs MC rapidity
+ TH2F* fph_reco_pMC_vs_yMC = nullptr; ///< Reconstructed track phase space (MC mom vs rapidity)
+ TH1F* fph_reco_purity = nullptr; ///< Purity of reconstructed tracks (\note purity requires MC information)
+
// Residuals and pulls at the first track point
TH1F* fph_fst_res_x = nullptr; ///< Residual of x at first track point [cm]
TH1F* fph_fst_res_y = nullptr; ///< Residual of y at first track point [cm]
diff --git a/reco/L1/qa/CbmCaTrackFitQa.cxx b/reco/L1/qa/CbmCaTrackFitQa.cxx
new file mode 100644
index 0000000000..7769050a36
--- /dev/null
+++ b/reco/L1/qa/CbmCaTrackFitQa.cxx
@@ -0,0 +1,122 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// @file CbmCaTrackFitQa.cxx
+/// @brief QA submodule for track fit results (residuals and puls) at selected z-coordinate (implementation)
+/// @since 03.04.2023
+/// @author Sergei Zharko <s.zharko@gsi.de>
+
+#include "CbmCaTrackFitQa.h"
+
+#include "CbmL1Track.h"
+
+#include "CaToolsMCData.h"
+#include "L1Field.h"
+#include "L1Fit.h"
+
+using cbm::ca::TrackFitQa;
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+TrackFitQa::TrackFitQa(const char* pointTag, const char* prefixName, TFolder* pParentFolder)
+ : CbmQaIO(pointTag, Form("%s_%s", prefixName, pointTag), pParentFolder)
+{
+ fStoringMode = EStoringMode::kSAMEDIR;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void TrackFitQa::Init()
+{
+ fph_res_x = MakeHisto<TH1F>("res_x", "", fBinsRESX, fLoRESX, fUpRESX);
+ fph_res_y = MakeHisto<TH1F>("res_y", "", fBinsRESY, fLoRESY, fUpRESY);
+ fph_res_t = MakeHisto<TH1F>("res_t", "", fBinsREST, fLoREST, fUpREST);
+ fph_res_tx = MakeHisto<TH1F>("res_tx", "", fBinsRESTX, fLoRESTX, fUpRESTX);
+ fph_res_ty = MakeHisto<TH1F>("res_ty", "", fBinsRESTY, fLoRESTY, fUpRESTY);
+ fph_res_qp = MakeHisto<TH1F>("res_qp", "", fBinsRESQP, fLoRESQP, fUpRESQP);
+ fph_res_vi = MakeHisto<TH1F>("res_vi", "", fBinsRESVI, fLoRESVI, fUpRESVI);
+
+ fph_pull_x = MakeHisto<TH1F>("pull_x", "", fBinsPULLX, fLoPULLX, fUpPULLX);
+ fph_pull_y = MakeHisto<TH1F>("pull_y", "", fBinsPULLY, fLoPULLY, fUpPULLY);
+ fph_pull_t = MakeHisto<TH1F>("pull_t", "", fBinsPULLT, fLoPULLT, fUpPULLT);
+ fph_pull_tx = MakeHisto<TH1F>("pull_tx", "", fBinsPULLTX, fLoPULLTX, fUpPULLTX);
+ fph_pull_ty = MakeHisto<TH1F>("pull_ty", "", fBinsPULLTY, fLoPULLTY, fUpPULLTY);
+ fph_pull_qp = MakeHisto<TH1F>("pull_qp", "", fBinsPULLQP, fLoPULLQP, fUpPULLQP);
+ fph_pull_vi = MakeHisto<TH1F>("pull_vi", "", fBinsPULLVI, fLoPULLVI, fUpPULLVI);
+
+ // Set histogram titles
+ TString sPrefix = (fsTitle.Length() > 0) ? fsTitle + " point " : "";
+ fph_res_x->SetTitle(sPrefix + " residual for x-coordinate;x_{reco} - x_{MC} [cm]");
+ fph_res_y->SetTitle(sPrefix + " residual for y-coordinate;y_{reco} - y_{MC} [cm]");
+ fph_res_t->SetTitle(sPrefix + " residual for time;t_{reco} - t_{MC} [ns]");
+ fph_res_tx->SetTitle(sPrefix + " residual for slope along x-axis;T_{x}^{reco} - T_{x}^{MC}");
+ fph_res_ty->SetTitle(sPrefix + " residual for slope along y-axis;T_{y}^{reco} - T_{y}^{MC}");
+ fph_res_qp->SetTitle(sPrefix + " residual for q/p;(q/p)_{reco} - (q/p)_{MC} [ec/GeV]");
+ fph_res_vi->SetTitle(sPrefix + " residual for inverse speed;v^{-1}_{reco} - v^{-1}_{MC} [c^{-1}]");
+
+ fph_pull_x->SetTitle(sPrefix + " pull for x-coordinate;(x_{reco} - x_{MC}) / #sigma_{x}");
+ fph_pull_y->SetTitle(sPrefix + " pull for y-coordinate;(y_{reco} - y_{MC}) / #sigma_{y}");
+ fph_pull_t->SetTitle(sPrefix + " pull for time;(t_{reco} - t_{MC}) / #sigma_{t}");
+ fph_pull_tx->SetTitle(sPrefix + " pull for slope along x-axis;(T_{x}^{reco} - T_{x}^{MC}) / #sigma_{T_{x}}");
+ fph_pull_ty->SetTitle(sPrefix + " pull for slope along y-axis;(T_{y}^{reco} - T_{y}^{MC}) / #sigma_{T_{y}}");
+ fph_pull_qp->SetTitle(sPrefix + " pull for q/p;((q/p)_{reco} - (q/p)_{MC}) / #sigma_{q/p}");
+ fph_pull_vi->SetTitle(sPrefix + " pull for inverse speed;(v^{-1}_{reco} - v^{-1}_{MC}) / #sigma_{v^{-1}}");
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void TrackFitQa::Fill(const L1TrackPar& trPar, const ::ca::tools::MCPoint& mcPoint, bool bTimeMeasured, double weight)
+{
+ // Probably, a bottleneck
+ L1Fit fitter;
+ fitter.SetParticleMass(fMass);
+ fitter.SetMask(fmask::One());
+ fitter.SetDoFitVelocity(true);
+ fitter.SetTrack(trPar);
+ L1FieldRegion fieldRegion;
+ fieldRegion.SetUseOriginalField(); // Precised magnetic field is used
+ fitter.Extrapolate(mcPoint.GetZOut(), fieldRegion);
+
+ const L1TrackPar& trParExtr = fitter.Tr(); // Track parameters extrapolated to given MC point
+
+ // ** Time-independent measurements **
+
+ double resX = trParExtr.GetX() - mcPoint.GetXOut(); // residual of x-position [cm]
+ double resY = trParExtr.GetY() - mcPoint.GetYOut(); // residual of y-position [cm]
+ double resTx = trParExtr.GetTx() - mcPoint.GetTxOut(); // residual of slope along x-axis [1]
+ double resTy = trParExtr.GetTy() - mcPoint.GetTyOut(); // residual of slope along y-axis [1]
+ double resQp = trParExtr.GetQp() - mcPoint.GetQpOut(); // residual of q/p [ec/GeV]
+
+ double pullX = resX / trParExtr.GetXErr(); // pull of x-position
+ double pullY = resY / trParExtr.GetYErr(); // pull of y-position
+ double pullTx = resTx / trParExtr.GetTxErr(); // pull of slope along x-axis
+ double pullTy = resTy / trParExtr.GetTyErr(); // pull of slope along y-axis
+ double pullQp = resQp / trParExtr.GetQpErr(); // pull of q/p
+
+ fph_res_x->Fill(resX, weight);
+ fph_res_y->Fill(resY, weight);
+ fph_res_tx->Fill(resTx, weight);
+ fph_res_ty->Fill(resTy, weight);
+ fph_res_qp->Fill(resQp, weight);
+
+ fph_pull_x->Fill(pullX, weight);
+ fph_pull_y->Fill(pullY, weight);
+ fph_pull_tx->Fill(pullTx, weight);
+ fph_pull_ty->Fill(pullTy, weight);
+ fph_pull_qp->Fill(pullQp, weight);
+
+
+ // ** Time-dependent measurements **
+ if (bTimeMeasured) {
+ double resT = trParExtr.GetTime() - mcPoint.GetTime(); // residual of time [ns]
+ double resVi = (trParExtr.GetInvSpeed() - mcPoint.GetInvSpeedOut()) * L1Constants::phys::kSpeedOfLightInv;
+ double pullT = resT / trParExtr.GetTimeErr(); // pull of time
+ double pullVi = resVi / trParExtr.GetInvSpeedErr(); // pull of inverse speed
+
+ fph_res_t->Fill(resT, weight);
+ fph_res_vi->Fill(resVi, weight);
+ fph_pull_t->Fill(pullT, weight);
+ fph_pull_vi->Fill(pullVi, weight);
+ }
+}
diff --git a/reco/L1/qa/CbmCaTrackFitQa.h b/reco/L1/qa/CbmCaTrackFitQa.h
new file mode 100644
index 0000000000..4f2ce84f43
--- /dev/null
+++ b/reco/L1/qa/CbmCaTrackFitQa.h
@@ -0,0 +1,163 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// @file CbmCaTrackFitQa.h
+/// @brief QA submodule for track fit results (residuals and puls) at selected z-coordinate (header)
+/// @since 03.04.2023
+/// @author Sergei Zharko <s.zharko@gsi.de>
+
+#ifndef CbmCaTrackFitQa_h
+#define CbmCaTrackFitQa_h 1
+
+#include "CbmQaIO.h"
+
+#include "L1Constants.h"
+#include "L1Field.h"
+#include "L1Fit.h"
+#include "L1Vector.h"
+
+// Forward declarations
+namespace ca::tools
+{
+ class MCData;
+ class MCPoint;
+} // namespace ca::tools
+class CbmL1Track;
+class CbmL1HitDebugInfo;
+
+class TFolder;
+class TH1F;
+
+namespace cbm::ca
+{
+ /// @brief Set of histograms to monitor track parameters
+ ///
+ /// Class provides histograms to monitor track fit parameters at a selected z-coordinate.
+ /// The parameters include x, y, tx, ty, time, q/p, vi (inverse speed).
+ class TrackFitQa : public CbmQaIO {
+ public:
+ /// @brief Constructor
+ /// @param pointTag Tag for point, in which the parameters are analyzed
+ /// @param prefixName Name of unique prefix
+ /// @param pParentFolder Pointer to parent folder
+ TrackFitQa(const char* pointTag, const char* prefixName, TFolder* pParentFolder);
+
+ /// @brief Destructor
+ ~TrackFitQa() = default;
+
+ /// @brief Copy constructor
+ TrackFitQa(const TrackFitQa&) = delete;
+
+ /// @brief Move constructor
+ TrackFitQa(TrackFitQa&&) = delete;
+
+ /// @brief Copy assignment operator
+ TrackFitQa& operator=(const TrackFitQa&) = delete;
+
+ /// @brief Move assignment operator
+ TrackFitQa& operator=(TrackFitQa&&) = delete;
+
+ /// @brief Gets title of fit parameters
+ const char* GetTitle() const { return fsTitle; }
+
+ /// @brief Initializes histograms
+ void Init();
+
+ /// @brief Fills pull and residual histograms
+ /// @note To be called in the loop over reconstructed tracks full sample
+ /// @param iTrkReco Index of reconstructed track
+ /// @param weight Weight
+ void Fill(const L1TrackPar& trPar, const ::ca::tools::MCPoint& mcPoint, bool bTimeMeasured, double weight = 1);
+
+ /// @brief Sets particle mass, used for fitting a track
+ /// @param mass Particle mass [GeV/c2]
+ void SetParticleMass(double mass) { fMass = mass; }
+
+ /// @brief Sets title, which is to be reflected on legends and titles
+ /// @param title Title of fit distributions
+ void SetTitle(const char* title) { fsTitle = title; }
+
+ // ************************
+ // ** List of histograms **
+ // ************************
+
+ // ** Residuals **
+ TH1F* fph_res_x = nullptr; ///< Residual of x-coordinate [cm]
+ TH1F* fph_res_y = nullptr; ///< Residual of y-coordinate [cm]
+ TH1F* fph_res_tx = nullptr; ///< Residual of slope along x-axis
+ TH1F* fph_res_ty = nullptr; ///< Residual of slope along y-axis
+ TH1F* fph_res_qp = nullptr; ///< Residual of q/p [ec/GeV]
+ TH1F* fph_res_t = nullptr; ///< Residual of time [ns]
+ TH1F* fph_res_vi = nullptr; ///< Residual of inverse speed [1/c]
+
+ // ** Pulls **
+ TH1F* fph_pull_x = nullptr; ///< Pull of x-coordinate
+ TH1F* fph_pull_y = nullptr; ///< Pull of y-coordinate
+ TH1F* fph_pull_tx = nullptr; ///< Pull of slope along x-axis
+ TH1F* fph_pull_ty = nullptr; ///< Pull of slope along y-axis
+ TH1F* fph_pull_qp = nullptr; ///< Pull of q/p
+ TH1F* fph_pull_t = nullptr; ///< Pull of time
+ TH1F* fph_pull_vi = nullptr; ///< Pull of inverse speed
+
+
+ // **************************
+ // ** Histogram properties **
+ // **************************
+
+ // ** Binning **
+ int fBinsRESX = 200; ///< Number of bins, residual of x
+ double fLoRESX = -4.e-3; ///< Lower boundary, residual of x [cm]
+ double fUpRESX = +4.e-3; ///< Upper boundary, residual of x [cm]
+ int fBinsRESY = 200; ///< Number of bins, residual of y
+ double fLoRESY = -4.e-2; ///< Lower boundary, residual of y [cm]
+ double fUpRESY = +4.e-2; ///< Upper boundary, residual of y [cm]
+ int fBinsRESTX = 200; ///< Number of bins, residual of slope along x-axis
+ double fLoRESTX = -4.e-3; ///< Lower boundary, residual of slope along x-axis
+ double fUpRESTX = +4.e-3; ///< Upper boundary, residual of slope along x-axis
+ int fBinsRESTY = 200; ///< Number of bins, residual of slope along y-axis
+ double fLoRESTY = -4.e-3; ///< Lower boundary, residual of slope along y-axis
+ double fUpRESTY = +4.e-3; ///< Upper boundary, residual of slope along y-axis
+ int fBinsRESQP = 200; ///< Number of bins, residual of q/p
+ double fLoRESQP = -10.; ///< Lower boundary, residual of q/p [ec/GeV]
+ double fUpRESQP = +10.; ///< Upper boundary, residual of q/p [ec/GeV]
+ int fBinsREST = 200; ///< Number of bins, residual of time
+ double fLoREST = -10.; ///< Lower boundary, residual of time [ns]
+ double fUpREST = +10.; ///< Upper boundary, residual of time [ns]
+ int fBinsRESVI = 200; ///< Number of bins, residual of inverse speed
+ double fLoRESVI = -2.; ///< Lower boundary, residual of inverse speed [1/c]
+ double fUpRESVI = +2.; ///< Upper boundary, residual of inverse speed [1/c]
+
+ int fBinsPULLX = 200; ///< Number of bins, pull of x
+ double fLoPULLX = -1.; ///< Lower boundary, pull of x [cm]
+ double fUpPULLX = +1.; ///< Upper boundary, pull of x [cm]
+ int fBinsPULLY = 200; ///< Number of bins, pull of y
+ double fLoPULLY = -1.; ///< Lower boundary, pull of y [cm]
+ double fUpPULLY = +1.; ///< Upper boundary, pull of y [cm]
+ int fBinsPULLTX = 200; ///< Number of bins, pull of slope along x-axis
+ double fLoPULLTX = -4.; ///< Lower boundary, pull of slope along x-axis
+ double fUpPULLTX = +4.; ///< Upper boundary, pull of slope along x-axis
+ int fBinsPULLTY = 200; ///< Number of bins, pull of slope along y-axis
+ double fLoPULLTY = -4.; ///< Lower boundary, pull of slope along y-axis
+ double fUpPULLTY = +4.; ///< Upper boundary, pull of slope along y-axis
+ int fBinsPULLQP = 200; ///< Number of bins, pull of q/p
+ double fLoPULLQP = -10.; ///< Lower boundary, pull of q/p [ec/GeV]
+ double fUpPULLQP = +10.; ///< Upper boundary, pull of q/p [ec/GeV]
+ int fBinsPULLT = 200; ///< Number of bins, pull of time
+ double fLoPULLT = -10.; ///< Lower boundary, pull of time [ns]
+ double fUpPULLT = +10.; ///< Upper boundary, pull of time [ns]
+ int fBinsPULLVI = 200; ///< Number of bins, pull of inverse speed
+ double fLoPULLVI = -2.; ///< Lower boundary, pull of inverse speed [1/c]
+ double fUpPULLVI = +2.; ///< Upper boundary, pull of inverse speed [1/c]
+
+ private:
+ TString fsTitle = ""; ///< Title of the point
+
+ double fMass = L1Constants::phys::kMuonMass; /// Mass of particle
+ };
+
+
+} // namespace cbm::ca
+
+
+#endif // CbmCaTrackFitQa_h
diff --git a/reco/L1/qa/CbmCaTrackTypeQa.cxx b/reco/L1/qa/CbmCaTrackTypeQa.cxx
new file mode 100644
index 0000000000..c69ce7aa16
--- /dev/null
+++ b/reco/L1/qa/CbmCaTrackTypeQa.cxx
@@ -0,0 +1,157 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// @file CbmCaTrackTypeQa.h
+/// @brief QA submodule for different track types (header)
+/// @since 27.03.2023
+/// @author Sergei Zharko <s.zharko@gsi.de>
+
+#include "CbmCaTrackTypeQa.h"
+
+#include "CbmCaTrackFitQa.h"
+#include "CbmL1Track.h"
+
+#include "CaToolsMCData.h"
+
+using ca::tools::MCPoint;
+using ca::tools::MCTrack;
+using cbm::ca::TrackTypeQa;
+
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+TrackTypeQa::TrackTypeQa(const char* typeName, const char* prefixName, bool bUseMC, TFolder* pParentFolder)
+ : CbmQaIO(typeName, Form("%s_%s", prefixName, typeName), pParentFolder)
+ , fbUseMC(bUseMC)
+{
+ fStoringMode = EStoringMode::kSAMEDIR;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void TrackTypeQa::Init()
+{
+ // TODO: Replace assertions with exceptions
+ assert(fpvRecoTracks);
+ assert(fpvHits);
+ assert(!fbUseMC || fpMCData);
+
+ fph_reco_p = MakeHisto<TH1F>("reco_p", "", kBinsP, kLoP, kUpP);
+ fph_reco_pt = MakeHisto<TH1F>("reco_pt", "", kBinsPT, kLoPT, kUpPT);
+ fph_reco_eta = MakeHisto<TH1F>("reco_eta", "", kBinsETA, kLoETA, kUpETA);
+ fph_reco_phi = MakeHisto<TH1F>("reco_phi", "", kBinsPHI, kLoPHI, kUpPHI);
+ fph_reco_theta = MakeHisto<TH1F>("reco_theta", "", kBinsTHETA, kLoTHETA, kUpTHETA);
+ fph_reco_tx = MakeHisto<TH1F>("reco_tx", "", kBinsTX, kLoTX, kUpTX);
+ fph_reco_ty = MakeHisto<TH1F>("reco_ty", "", kBinsTY, kLoTY, kUpTY);
+ fph_reco_fhitR = MakeHisto<TH1F>("reco_fhitR", "", kBinsFHITR, kLoFHITR, kUpFHITR);
+ fph_reco_nhits = MakeHisto<TH1F>("reco_nhits", "", kBinsNHITS, kLoNHITS, kUpNHITS);
+ // TODO: ...
+
+ fph_reco_p->SetTitle("Total momentum of reconstructed track;p_{reco} [GeV/c]");
+ fph_reco_pt->SetTitle("Transverse momentum of reconstructed track;p_{T}_{reco} [GeV/c]");
+ fph_reco_phi->SetTitle("Azimuthal angle of reconstructed track;#phi_{reco} [rad]");
+ fph_reco_theta->SetTitle("Polar angle of reconstructed track;#theta_{reco} [rad]");
+ fph_reco_tx->SetTitle("Slope along x-axis of reconstructed tracks;t_{x}");
+ fph_reco_ty->SetTitle("Slope along y-axis of reconstructed tracks;t_{y}");
+ fph_reco_fhitR->SetTitle("Distance of the first hit from z-axis for reconstructed tracks");
+ fph_reco_nhits->SetTitle("Hit number of reconstructed tracks");
+
+ if (fbUseMC) {
+ fph_reco_pMC = MakeHisto<TH1F>("reco_pMC", "", kBinsP, kLoP, kUpP);
+ fph_reco_yMC = MakeHisto<TH1F>("reco_yMC", "", kBinsY, kLoY, kUpY);
+ fph_reco_pMC_yMC = MakeHisto<TH2F>("reco_pMC_yMC", "", kBinsY, kLoY, kUpY, kBinsP, kLoP, kUpP);
+ fph_reco_phiMC = MakeHisto<TH1F>("reco_phiMC", "", kBinsPHI, kLoPHI, kUpPHI);
+ fph_reco_thetaMC = MakeHisto<TH1F>("reco_thetaMC", "", kBinsTHETA, kLoTHETA, kUpTHETA);
+
+ fph_reco_pMC->SetTitle("MC total momentum of reconstructed track;p_{MC} [GeV/c]");
+ fph_reco_yMC->SetTitle("MC rapidity of reconstructed track;y_{MC}");
+ fph_reco_pMC_yMC->SetTitle("Transverse momentum of reconstructed track;y_{MC};p_{T}_{MC} [GeV/c]");
+ fph_reco_phiMC->SetTitle("Azimuthal angle of reconstructed track;#phi_{reco} [rad]");
+ fph_reco_thetaMC->SetTitle("Polar angle of reconstructed track;#theta_{reco} [rad]");
+
+ fph_mc_pMC = MakeHisto<TH1F>("mc_pMC", "", kBinsP, kLoP, kUpP);
+ fph_mc_yMC = MakeHisto<TH1F>("mc_yMC", "", kBinsY, kLoY, kUpY);
+ fph_mc_pMC_yMC = MakeHisto<TH2F>("mc_pMC_yMC", "", kBinsY, kLoY, kUpY, kBinsP, kLoP, kUpP);
+
+ fph_mc_pMC->SetTitle("MC total momentum of MC tracks;p_{MC} [GeV/c]");
+ fph_mc_yMC->SetTitle("MC rapidity of MC tracks;y_{MC}");
+ fph_mc_pMC_yMC->SetTitle("MC total momentum vs. MC rapidity of MC tracks;y_{MC};p_{MC} [GeV/c]");
+
+ fpFitQaFirstHit = std::make_unique<TrackFitQa>("fst_hit", fsPrefix, fpFolderRoot);
+ fpFitQaFirstHit->SetTitle("First hit");
+ fpFitQaFirstHit->Init();
+
+ fpFitQaLastHit = std::make_unique<TrackFitQa>("lst_hit", fsPrefix, fpFolderRoot);
+ fpFitQaLastHit->SetTitle("Last hit");
+ fpFitQaLastHit->Init();
+
+ fpFitQaVertex = std::make_unique<TrackFitQa>("vertex", fsPrefix, fpFolderRoot);
+ fpFitQaVertex->SetTitle("Vertex");
+ fpFitQaVertex->Init();
+ }
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void TrackTypeQa::FillRecoTrack(int iTrkReco, double weight)
+{
+ const auto& recoTrack = (*fpvRecoTracks)[iTrkReco];
+ fph_reco_p->Fill(recoTrack.GetP(), weight);
+ fph_reco_pt->Fill(recoTrack.GetPt(), weight);
+ fph_reco_phi->Fill(recoTrack.GetPhi(), weight);
+ fph_reco_theta->Fill(recoTrack.GetTheta(), weight);
+ fph_reco_tx->Fill(recoTrack.GetTx(), weight);
+ fph_reco_ty->Fill(recoTrack.GetTy(), weight);
+ if (fbUseMC) {
+ int iTrkMC = recoTrack.GetMatchedMCTrackIndex();
+ if (iTrkMC > -1) {
+ /// TODO: fill mass hypothesis into CbmL1Track
+ const auto& mcTrack = fpMCData->GetTrack(iTrkMC);
+ fph_reco_pMC->Fill(mcTrack.GetP(), weight);
+ fph_reco_yMC->Fill(mcTrack.GetRapidity(), weight);
+ fph_reco_pMC_yMC->Fill(mcTrack.GetRapidity(), mcTrack.GetP(), weight);
+
+ // *****************************
+ // ** Track fit parameters QA **
+ // *****************************
+
+ int nTimeMeasurements = 0;
+ for (int iH : recoTrack.GetHitIndexes()) {
+ int iSt = (*fpvHits)[iH].GetStationId();
+ if (fpParameters->GetStation(iSt).timeInfo) { ++nTimeMeasurements; }
+ }
+ bool isTimeFitted = (nTimeMeasurements > 1);
+
+ // ** First hit **
+ int iHfst = recoTrack.GetFirstHitIndex();
+ int iPfst = (*fpvHits)[iHfst].GetMCPointIndex();
+ if (iPfst > -1) {
+ const auto& mcPoint = fpMCData->GetPoint(iPfst);
+ L1TrackPar trPar(recoTrack.T, recoTrack.C);
+ fpFitQaFirstHit->Fill(trPar, mcPoint, isTimeFitted, weight);
+ }
+
+ // ** Last hit **
+ int iHlst = recoTrack.GetLastHitIndex();
+ int iPlst = (*fpvHits)[iHlst].GetMCPointIndex();
+ if (iPlst > -1) {
+ const auto& mcPoint = fpMCData->GetPoint(iPlst);
+ L1TrackPar trPar(recoTrack.TLast, recoTrack.CLast);
+ fpFitQaLastHit->Fill(trPar, mcPoint, isTimeFitted, weight);
+ }
+ }
+ }
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void TrackTypeQa::FillMCTrack(int iTrkMC, double weight)
+{
+ assert(fbUseMC);
+
+ const MCTrack& mcTrack = fpMCData->GetTrack(iTrkMC);
+ fph_mc_pMC->Fill(mcTrack.GetP(), weight);
+ fph_mc_yMC->Fill(mcTrack.GetRapidity(), weight);
+ fph_mc_pMC_yMC->Fill(mcTrack.GetRapidity(), mcTrack.GetP(), weight);
+}
diff --git a/reco/L1/qa/CbmCaTrackTypeQa.h b/reco/L1/qa/CbmCaTrackTypeQa.h
new file mode 100644
index 0000000000..c54479019f
--- /dev/null
+++ b/reco/L1/qa/CbmCaTrackTypeQa.h
@@ -0,0 +1,209 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// @file CbmCaTrackTypeQa.h
+/// @brief QA submodule for different track types (header)
+/// @since 27.03.2023
+/// @author Sergei Zharko <s.zharko@gsi.de>
+
+#ifndef CbmCaTrackTypeQa_h
+#define CbmCaTrackTypeQa_h 1
+
+#include "CbmCaTrackFitQa.h"
+#include "CbmL1Hit.h"
+#include "CbmQaIO.h"
+
+#include <map>
+#include <string>
+
+#include "L1Parameters.h"
+
+// Forward declarations
+namespace ca::tools
+{
+ class MCData;
+}
+class CbmL1Track;
+class CbmL1HitDebugInfo;
+
+class TH1F;
+class TH2F;
+class TProfile;
+
+namespace cbm::ca
+{
+ /// @brief Unified QA for a group of tracks
+ ///
+ /// The class provides set of histograms, counters, and efficiencies for monitoring tracks of a particular type
+ /// (reco_fast_long, mc_pi, reco_ghost)
+ class TrackTypeQa : public CbmQaIO {
+ public:
+ /// @brief Constructor
+ /// @param typeName Name of tracks type (used as a prefix for histogram)
+ /// @note Track type is stored as a root directory name in the CbmQaIO base class
+ /// @param prefixName Name of unique prefix
+ /// @param bUseMC Flag: true - MC is used
+ /// @param pParentFolder Pointer to parent folder
+ TrackTypeQa(const char* typeName, const char* prefixName, bool bUseMC, TFolder* pParentFolder);
+
+ /// @brief Destructor
+ ~TrackTypeQa() = default;
+
+ /// @brief Copy constructor
+ TrackTypeQa(const TrackTypeQa&) = delete;
+
+ /// @brief Move constructor
+ TrackTypeQa(TrackTypeQa&&) = delete;
+
+ /// @brief Copy assignment operator
+ TrackTypeQa& operator=(const TrackTypeQa&) = delete;
+
+ /// @brief Move assignment operator
+ TrackTypeQa& operator=(TrackTypeQa&&) = delete;
+
+ /// @brief Gets title of track group
+ const char* GetTitle() const { return fsTitle; }
+
+ /// @brief Initializes histograms
+ void Init();
+
+ /// @brief Fills histograms with various track information
+ /// @note To be called in the loop over reconstructed tracks full sample
+ /// @param iTrkReco Index of reconstructed track
+ /// @param weight Weight
+ void FillRecoTrack(int iTrkReco, double weight = 1);
+
+ /// @brief Fills histograms with mc track information
+ /// @note To be called in the loop over MC tracks full sample
+ /// @param iTrkMC Index of MC track
+ /// @param weight Weight
+ void FillMCTrack(int iTrkMC, double weight = 1);
+
+ /// @brief Registers the reconstructed tracks container
+ /// @param vTracks Reference to reconstructed tracks container
+ void RegisterRecoTracks(L1Vector<CbmL1Track>& vTracks) { fpvRecoTracks = &vTracks; }
+
+ /// @brief Register the reconstructed hits container
+ /// @param vHits Reference to reconstructed hits container
+ void RegisterRecoHits(L1Vector<CbmL1HitDebugInfo>& vHits) { fpvHits = &vHits; }
+
+ /// @brief Registers the MC data
+ /// @param mcData Reference to MC data object
+ void RegisterMCData(::ca::tools::MCData& mcData) { fpMCData = &mcData; }
+
+ /// @brief Registers CA parameters object
+ /// @param pParameters A shared pointer to the parameters object
+ void RegisterParameters(std::shared_ptr<L1Parameters>& pParameters) { fpParameters = pParameters; }
+
+ /// @brief Sets title, which is to be reflected on legends and titles
+ /// @param title Title of track type
+ void SetTitle(const char* title) { fsTitle = title; }
+
+ // ************************
+ // ** List of histograms **
+ // ************************
+
+ // NOTE: Naming convention:
+ // "fph_" + <"reco"/"mc"> + <quantity> + <""/"MC">. Here <"reco"/"mc"> stands for reconstructed track or true MC
+ // track, <quantity> is a name of a quantity, versus which one wants to build the dependency, <""/"MC"> stands
+ // for an origin of the quantity - reconstructed or true Monte-Carlo. If there are several quantities, they are
+ // to be separated with "_" (example: "p_yMC" -- reconstructed momentum vs MC rapidity)
+
+ // ** Histograms from reconstructed information only **
+ TH1F* fph_reco_p = nullptr; ///< Total momentum over charge of reconstructed tracks
+ TH1F* fph_reco_pt = nullptr; ///< Transverse momentum over charge of reconstructed tracks
+ TH1F* fph_reco_phi = nullptr; ///< Azimuthal angle of reconstructed tracks
+ TH1F* fph_reco_theta = nullptr; ///< Polar angle of reconstructed tracks
+ TH1F* fph_reco_tx = nullptr; ///< Slope along x-axis of reconstructed tracks
+ TH1F* fph_reco_ty = nullptr; ///< Slope along y-axis of reconstructed tracks
+ TH1F* fph_reco_eta = nullptr; ///< Pseudo-rapidity of reconstructed tracks
+ TH1F* fph_reco_fhitR = nullptr; ///< Distance of the first hit from z-axis for reconstructed tracks
+ TH1F* fph_reco_nhits = nullptr; ///< Hit number of reconstructed tracks
+ TH1F* fph_reco_fsta = nullptr; ///< First station index of reconstructed tracks
+ TH1F* fph_reco_chi2_ndf = nullptr; ///< Fit chi2 over NDF of reconstructed tracks
+
+ // ** Reconstructed track distributions, requiring MC information **
+ TH1F* fph_reco_pMC = nullptr; ///< MC total momentum over charge of reconstructed tracks
+ TH1F* fph_reco_yMC = nullptr; ///< MC rapidity of reconstructed tracks
+ TH1F* fph_reco_etaMC = nullptr; ///< MC pseudo-rapidity of reconstructed tracks
+ TH2F* fph_reco_p_yMC = nullptr; ///< Total momentum vs MC rapidity of reconstructed tracks
+ TH2F* fph_reco_pMC_yMC = nullptr; ///< MC total momentum vs MC rapidity of reconstructed tracks
+ TH1F* fph_reco_ptMC = nullptr; ///< Transverse momentum over charge of reconstructed tracks
+ TH1F* fph_reco_phiMC = nullptr; ///< Azimuthal angle of reconstructed tracks
+ TH1F* fph_reco_thetaMC = nullptr; ///< Polar angle of reconstructed tracks
+ TH1F* fph_reco_txMC = nullptr; ///< Slope along x-axis of reconstructed tracks
+ TH1F* fph_reco_tyMC = nullptr; ///< Slope along y-axis of reconstructed tracks
+
+ // ** MC track distributions **
+ TH1F* fph_mc_pMC = nullptr; ///< MC total momentum over charge of MC tracks
+ TH1F* fph_mc_yMC = nullptr; ///< MC rapidity of MC tracks
+ TH1F* fph_mc_etaMC = nullptr; ///< MC pseudo-rapidity of MC tracks
+ TH2F* fph_mc_pMC_yMC = nullptr; ///< MC total momentum vs. MC rapidity of MC tracks
+ TH1F* fph_mc_ptMC = nullptr; ///< Transverse momentum over charge of reconstructed tracks
+ TH1F* fph_mc_phiMC = nullptr; ///< Azimuthal angle of reconstructed tracks
+ TH1F* fph_mc_thetaMC = nullptr; ///< Polar angle of reconstructed tracks
+ TH1F* fph_mc_txMC = nullptr; ///< Slope along x-axis of reconstructed tracks
+ TH1F* fph_mc_tyMC = nullptr; ///< Slope along y-axis of reconstructed tracks
+
+ // ** Efficiencies **
+ TProfile* fph_eff_p = nullptr; ///< Track efficiency
+ TProfile* fph_eff_nhits = nullptr;
+ TProfile* fph_eff_pMC = nullptr;
+
+ // ** Fit QA **
+ std::unique_ptr<TrackFitQa> fpFitQaFirstMCpoint = nullptr;
+ std::unique_ptr<TrackFitQa> fpFitQaFirstHit = nullptr;
+ std::unique_ptr<TrackFitQa> fpFitQaLastHit = nullptr;
+ std::unique_ptr<TrackFitQa> fpFitQaVertex = nullptr;
+
+ // TODO: Provide integrated efficiencies
+
+ private:
+ bool fbUseMC = false; ///< Flag: true - MC information is used
+ TString fsTitle = ""; ///< Title of the track category
+
+ L1Vector<CbmL1Track>* fpvRecoTracks = nullptr; ///< Pointer to vector of reconstructed tracks
+ L1Vector<CbmL1HitDebugInfo>* fpvHits = nullptr; ///< Pointer to vector of reconstructed hits
+ ::ca::tools::MCData* fpMCData = nullptr; ///< Pointer to MC data object
+ std::shared_ptr<L1Parameters> fpParameters = nullptr; ///< Pointer to parameters object
+
+ // **************************
+ // ** Histogram properties **
+ // **************************
+
+ // ** Binning **
+ static constexpr int kBinsP = 150; ///< Number of bins, total momentum
+ static constexpr double kLoP = 0.; ///< Lower boundary, total momentum [GeV/c]
+ static constexpr double kUpP = 15.; ///< Upper boundary, total momentum [GeV/c]
+ static constexpr int kBinsPT = 100; ///< Number of bins, transverse momentum
+ static constexpr double kLoPT = 0.; ///< Lower boundary, transverse momentum [GeV/c]
+ static constexpr double kUpPT = 10.; ///< Upper boundary, transverse momentum [GeV/c]
+ static constexpr int kBinsETA = 40; ///< Number of bins, pseudo-rapidity
+ static constexpr double kLoETA = 0.; ///< Lower boundary, pseudo-rapidity
+ static constexpr double kUpETA = 4.; ///< Upper boundary, pseudo-rapidity
+ static constexpr int kBinsY = 40; ///< Number of bins, rapidity
+ static constexpr double kLoY = 0.; ///< Lower boundary, rapidity
+ static constexpr double kUpY = 4.; ///< Upper boundary, rapidity
+ static constexpr int kBinsPHI = 68; ///< Number of bins, azimuthal angle
+ static constexpr double kLoPHI = -3.2; ///< Lower boundary, azimuthal angle [rad]
+ static constexpr double kUpPHI = +3.2; ///< Upper boundary, azimuthal angle [rad]
+ static constexpr int kBinsTHETA = 68; ///< Number of bins, polar angle
+ static constexpr double kLoTHETA = 0.; ///< Lower boundary, polar angle [rad]
+ static constexpr double kUpTHETA = 3.2; ///< Upper boundary, polar angle [rad]
+ static constexpr int kBinsTX = 80; ///< Number of bins, slope along x
+ static constexpr double kLoTX = -2.; ///< Lower boundary, slope along x
+ static constexpr double kUpTX = +2.; ///< Upper boundary, slope along x
+ static constexpr int kBinsTY = 80; ///< Number of bins, slope along y
+ static constexpr double kLoTY = -2.; ///< Lower boundary, slope along y
+ static constexpr double kUpTY = +2.; ///< Upper boundary, slope along y
+ static constexpr int kBinsNHITS = 15; ///< Number of bins, number of hits
+ static constexpr double kLoNHITS = -0.5; ///< Lower boundary, number of hits
+ static constexpr double kUpNHITS = 14.5; ///< Upper boundary, number of hits
+ static constexpr int kBinsFHITR = 50; ///< Number of bins, transverse distance of the 1st hit from z-axis
+ static constexpr double kLoFHITR = 0.; ///< Lower boundary, transverse distance of the 1st hit from z-axis [cm]
+ static constexpr double kUpFHITR = 150.; ///< Upper boundary, transverse distance of the 1st hit from z-axis [cm]
+ };
+} // namespace cbm::ca
+
+#endif // CbmCaOutput
\ No newline at end of file
--
GitLab