diff --git a/core/qa/CMakeLists.txt b/core/qa/CMakeLists.txt
index f59b755abec2cd102a6f2a42eea4916fd354ffae..f71105c8a624183d6d4f72d61b61c7eb287ab708 100644
--- a/core/qa/CMakeLists.txt
+++ b/core/qa/CMakeLists.txt
@@ -17,6 +17,7 @@ set(LINKDEF ${LIBRARY_NAME}LinkDef.h)
set(PUBLIC_DEPENDENCIES
FairRoot::Base
FairLogger::FairLogger
+ external::yaml-cpp
ROOT::Core
ROOT::Gpad
ROOT::Hist
diff --git a/core/qa/CbmQaTask.cxx b/core/qa/CbmQaTask.cxx
index a534f1758b23431830a0b34e19d833585f5d3f82..aed9e84f6c4983f09a51c88bab5822657fa67f06 100644
--- a/core/qa/CbmQaTask.cxx
+++ b/core/qa/CbmQaTask.cxx
@@ -65,9 +65,24 @@ void CbmQaTask::Finish()
InitStatus CbmQaTask::Init()
{
LOG_IF(info, fVerbose > 0) << fName << ": initializing task ...";
-
InitStatus res = kSUCCESS;
+ // ----- Open configuration file
+ YAML::Node config;
+ if (fsConfigName.Length()) {
+ LOG(info) << fName << ": reading configuration from file \"" << fsConfigName << "\"";
+ try {
+ config = YAML::LoadFile(fsConfigName.Data())["qa"][fName.Data()];
+ fpCurrentNode = &config;
+ }
+ catch (const YAML::BadFile& exc) {
+ LOG(fatal) << fName << ": configuration file \"" << fsConfigName << "\" does not exist";
+ }
+ catch (const YAML::ParserException& exc) {
+ LOG(fatal) << fName << ": configuration file \"" << fsConfigName << "\" is formatted improperly";
+ }
+ }
+
// ----- Clear map of the histograms (note)
DeInitBase();
@@ -80,12 +95,13 @@ InitStatus CbmQaTask::Init()
fpFolderRoot->SetOwner(true); // When true, TFolder owns all added objects
fpFolderRoot->Add(&fNofEvents);
- // ----- Initialize histograms and canvases
+ // ----- Initialize histograms
LOG_IF(info, fVerbose > 1) << fName << ": initializing histograms";
res = std::max(res, InitHistograms());
fNofEvents.SetVal(0);
+ fpCurrentNode = nullptr; // De-init pointer to
return res;
}
@@ -111,57 +127,3 @@ void CbmQaTask::DeInitBase()
// De-initialize particular QA-task implementation
DeInit();
}
-
-
-// ---------------------------------------------------------------------------------------------------------------------
-// NOTE: Example
-void CbmQaTask::SetHistoProperties(TH1* pHisto)
-{
- constexpr double kHdefTextSize = 0.04;
-
- pHisto->SetStats(kTRUE);
- pHisto->Sumw2();
-
- // Axis property settings
- std::array<TAxis*, 3> vpAxes = {pHisto->GetXaxis(), pHisto->GetYaxis(), pHisto->GetZaxis()};
- for (auto* pAxis : vpAxes) {
- pAxis->SetTitleSize(kHdefTextSize);
- pAxis->SetLabelSize(kHdefTextSize);
- }
- vpAxes[0]->SetNdivisions(504);
-}
-
-// ---------------------------------------------------------------------------------------------------------------------
-//
-InitStatus CbmQaTask::InitDataBranches()
-{
- LOG_IF(info, fVerbose > 1) << fName << ": data branches initialization function is not defined";
- return kSUCCESS;
-}
-
-// ---------------------------------------------------------------------------------------------------------------------
-//
-InitStatus CbmQaTask::InitHistograms()
-{
- LOG_IF(info, fVerbose > 1) << fName << ": histogram initialization function is not defined";
- return kSUCCESS;
-}
-
-// ---------------------------------------------------------------------------------------------------------------------
-//
-void CbmQaTask::FillHistograms()
-{
- LOG_IF(info, fVerbose > 1) << fName << ": histogram filling function is not defined";
-}
-
-// ---------------------------------------------------------------------------------------------------------------------
-//
-InitStatus CbmQaTask::InitCanvases()
-{
- LOG_IF(info, fVerbose > 1) << fName << ": no initialization of canvases is provided";
- return kSUCCESS;
-}
-
-// ---------------------------------------------------------------------------------------------------------------------
-//
-void CbmQaTask::DeInit() { LOG_IF(info, fVerbose > 1) << fName << ": no extra de-initialization is provided"; }
diff --git a/core/qa/CbmQaTask.h b/core/qa/CbmQaTask.h
index 61f876d7e9af887b9f5c8c73b6f648868e566a1a..143b7a7abde728c1c6256a08c8371652f72a048d 100644
--- a/core/qa/CbmQaTask.h
+++ b/core/qa/CbmQaTask.h
@@ -19,11 +19,20 @@
#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 <algorithm>
+#include <regex>
#include <string_view>
+#include <type_traits>
+
+#include <yaml-cpp/yaml.h>
/// Class CbmQaTask is to be inherited with a particular QA-task. It provides mechanisms for storage and management
/// of QA canvases and histograms management
@@ -63,6 +72,14 @@ public:
/// FairTask: Defines action of the task in the end of run
void Finish();
+ /// @brief Sets task configuration file
+ /// @param filename Name of file
+ ///
+ /// Sets name of YAML configuration file, which defines parameters needed for the task. The file is assumed
+ /// to contain root branch "qa/[task name]". A subbranch "histograms" contains a list of histograms with pre defined
+ /// parameters: title, number of bins and axis ranges.
+ void SetConfigName(const char* filename) { fsConfigName = filename; }
+
ClassDef(CbmQaTask, 0);
protected:
@@ -74,22 +91,22 @@ protected:
virtual bool Check() = 0;
/// De-initialize the task
- virtual void DeInit();
+ virtual void DeInit() {}
/// Initializes data branches
- virtual InitStatus InitDataBranches();
+ virtual InitStatus InitDataBranches() { return kSUCCESS; }
/// Initializes histograms
- virtual InitStatus InitHistograms();
+ virtual InitStatus InitHistograms() { return kSUCCESS; }
/// Initializes canvases
- virtual InitStatus InitCanvases();
+ virtual InitStatus InitCanvases() { return kSUCCESS; }
/// Initializes event / time-slice
virtual InitStatus InitTimeSlice() { return kSUCCESS; }
/// Method to fill histograms per event or time-slice
- virtual void FillHistograms();
+ 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)
@@ -112,22 +129,22 @@ protected:
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
+ /// @param id0 First index (optional)
+ /// @param id1 Second index (optional)
+ /// @param id2 Third index (optional)
+ /// @return Pointer to the histogram object
+ 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);
- /// Method to setup histogram properties (text sizes etc.)
- /// \param pHist Pointer to a histogram to tune
- virtual void SetHistoProperties(TH1* pHist);
-
- /// Gets a reference to histograms map
- ///const auto& GetHistoMap() const { return fmpHistList; }
-
- /// Gets a reference to canvases map
- ///const auto& GetCanvasMap() const { return fmpCanvList; }
-
/// Get current event number
int GetEventNumber() const { return fNofEvents.GetVal(); }
@@ -147,17 +164,9 @@ protected:
private:
- // **********************************************************
- // ** Functions accessible inside the CbmQaTask class only **
- // **********************************************************
-
- /// De-initialize this task
+ /// @brief De-initializes this task
void DeInitBase();
- // ***********************
- // ** Private variables **
- // ***********************
-
// I/O
std::unique_ptr<TFolder> fpFolderRoot = nullptr; ///< Root folder to store histograms and canvases
TFolder* fpFolderHist = nullptr; ///< Folder for raw histograms
@@ -169,6 +178,9 @@ private:
bool fbUseMC = false; ///< Flag, if MC is used
TParameter<int> fNofEvents {"nEvents", 0}; ///< Number of processed events
+
+ TString fsConfigName = ""; ///< Name of YAML configuration file
+ YAML::Node* fpCurrentNode = nullptr; ///< Pointer to current node of the configuration file
};
@@ -256,14 +268,131 @@ T* CbmQaTask::MakeHisto(const char* nameBase, Args... args)
delete pHist;
}
- // Create a new histogram or profile
T* pHist = new T(name, args...);
// Register histogram in the folder
if (!fpFolderHist) { fpFolderHist = fpFolderRoot->AddFolder("histograms", "Histograms"); }
fpFolderHist->Add(pHist);
- SetHistoProperties(pHist);
+ return pHist;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+template<typename T>
+T* CbmQaTask::MakeHistoFromConfig(const char* nameBase, int id0, int id1, int id2)
+{
+ std::string name = std::string(fsPrefix.Data()) + "_" + nameBase;
+ if (id0 > -1) { name = std::regex_replace(name, std::regex("\\%0"), std::to_string(id0)); }
+ if (id1 > -1) { name = std::regex_replace(name, std::regex("\\%1"), std::to_string(id1)); }
+ if (id2 > -1) { name = std::regex_replace(name, std::regex("\\%2"), std::to_string(id2)); }
+
+ // Check, if the histogram with a given name was already created. If so, delete it
+ if (gROOT->FindObject(name.data())) {
+ 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.data());
+ delete pHist;
+ }
+
+ T* pHist = nullptr;
+ // Create histogram
+ LOG_IF(fatal, !fpCurrentNode)
+ << fName << ": attempt to make a histogram (\"" << nameBase << "\") from configuration file, which was not "
+ << "defined. Please, provide configuration file with defined histogram list to the task via "
+ << "SetConfigName(filename) function.";
+ try {
+ const auto& node = (*fpCurrentNode)["histograms"][nameBase];
+ LOG_IF(fatal, !node) << fName << ": node for histogram \"" << nameBase
+ << "\" was not defined in the configuration file";
+
+ std::string title = node["title"].as<std::string>("").data();
+ if (id0 > -1) { title = std::regex_replace(title, std::regex("\\%0"), std::to_string(id0)); }
+ if (id1 > -1) { title = std::regex_replace(title, std::regex("\\%1"), std::to_string(id1)); }
+ if (id2 > -1) { title = std::regex_replace(title, std::regex("\\%2"), std::to_string(id2)); }
+
+ // 1D-profiles
+ if constexpr (std::is_base_of_v<TProfile, T>) {
+ int nBinsX = node["nbinsx"].as<int>();
+ double minX = node["minx"].as<double>();
+ double maxX = node["maxx"].as<double>();
+ double minY = node["miny"].as<double>(0.);
+ double maxY = node["maxy"].as<double>(0.);
+ std::string opt = node["opt"].as<std::string>("");
+ pHist = new T(name.data(), title.data(), nBinsX, minX, maxX, minY, maxY, opt.data());
+ }
+ // 2D-profiles
+ else if constexpr (std::is_base_of_v<TProfile2D, T>) {
+ int nBinsX = node["nbinsx"].as<int>();
+ double minX = node["minx"].as<double>();
+ double maxX = node["maxx"].as<double>();
+ int nBinsY = node["nbinsy"].as<int>();
+ double minY = node["miny"].as<double>();
+ double maxY = node["maxy"].as<double>();
+ double minZ = node["minz"].as<double>(0.);
+ double maxZ = node["maxz"].as<double>(0.);
+ std::string opt = node["opt"].as<std::string>("");
+ pHist = new T(name.data(), title.data(), nBinsX, minX, maxX, nBinsY, minY, maxY, minZ, maxZ, opt.data());
+ }
+ // 3D-profiles
+ else if constexpr (std::is_base_of_v<TProfile3D, T>) {
+ int nBinsX = node["nbinsx"].as<int>();
+ double minX = node["minx"].as<double>();
+ double maxX = node["maxx"].as<double>();
+ int nBinsY = node["nbinsy"].as<int>();
+ double minY = node["miny"].as<double>();
+ double maxY = node["maxy"].as<double>();
+ int nBinsZ = node["nbinsz"].as<int>();
+ double minZ = node["minz"].as<double>();
+ double maxZ = node["maxz"].as<double>();
+ std::string opt = node["opt"].as<std::string>("");
+ pHist = new T(name.data(), title.data(), nBinsX, minX, maxX, nBinsY, minY, maxY, minZ, maxZ, opt.data());
+ }
+ // 2D-histograms
+ else if constexpr (std::is_base_of_v<TH2, T>) {
+ int nBinsX = node["nbinsx"].as<int>();
+ double minX = node["minx"].as<double>();
+ double maxX = node["maxx"].as<double>();
+ int nBinsY = node["nbinsy"].as<int>();
+ double minY = node["miny"].as<double>();
+ double maxY = node["maxy"].as<double>();
+ pHist = new T(name.data(), title.data(), nBinsX, minX, maxX, nBinsY, minY, maxY);
+ }
+ // 3D-histograms + derived
+ else if constexpr (std::is_base_of_v<TH3, T>) {
+ int nBinsX = node["nbinsx"].as<int>();
+ double minX = node["minx"].as<double>();
+ double maxX = node["maxx"].as<double>();
+ int nBinsY = node["nbinsy"].as<int>();
+ double minY = node["miny"].as<double>();
+ double maxY = node["maxy"].as<double>();
+ int nBinsZ = node["nbinsz"].as<int>();
+ double minZ = node["minz"].as<double>();
+ double maxZ = node["maxz"].as<double>();
+ pHist = new T(name.data(), title.data(), nBinsX, minX, maxX, nBinsY, minY, maxY, nBinsZ, minZ, maxZ);
+ }
+ // 1D-histograms + derived
+ else if constexpr (std::is_base_of_v<TH1, T>) {
+ int nBinsX = node["nbinsx"].as<int>();
+ double minX = node["minx"].as<double>();
+ double maxX = node["maxx"].as<double>();
+ pHist = new T(name.data(), title.data(), nBinsX, minX, maxX);
+ }
+ // Other types are restricted
+ //else {
+ // static_assert(false, "CbmQaTask::MakeTable: unsupported type given as a template parameter");
+ //}
+ }
+ catch (const YAML::InvalidNode& exc) {
+ LOG(fatal) << fName << ": error while reading the histogram \"" << nameBase << "\" properties from "
+ << "configuration file \"" << fsConfigName << "\". Please, check the file formatting. \n Tip: probably, "
+ << "there ara mistakes in the obligatory property names, for example, \"nbinsX\" or \"nXBins\" is "
+ << "used instead of \"nbinsx\" etc.";
+ }
+
+ // Register histogram in the folder
+ if (!fpFolderHist) { fpFolderHist = fpFolderRoot->AddFolder("histograms", "Histograms"); }
+ fpFolderHist->Add(pHist);
return pHist;
}
diff --git a/macro/run/CMakeLists.txt b/macro/run/CMakeLists.txt
index f2965419df60ef7f15b59fe154155c09a9582271..c8846808a37c387d89f6f59384823eee4c792f82 100644
--- a/macro/run/CMakeLists.txt
+++ b/macro/run/CMakeLists.txt
@@ -15,6 +15,11 @@ file(COPY ${CBMROOT_SOURCE_DIR}/macro/include/.rootrc
DESTINATION ${CBMROOT_BINARY_DIR}/macro/run)
# ============================================================================
+# ===== Copy configuration file for QA tasks to the build directory ========
+file(COPY ${CBMROOT_SOURCE_DIR}/macro/run/qa_config.cbm.yaml
+ DESTINATION ${CBMROOT_BINARY_DIR}/macro/run)
+# ============================================================================
+
# ===== Define variables for tests =======================================
if(${CBM_TEST_MODEL} MATCHES Weekly)
@@ -223,7 +228,7 @@ foreach(setup IN LISTS cbm_setup)
add_test(${testname} ${MACRODIR}/run_qa.sh
\"data/${sname}_coll\" \"data/${sname}_ts\" \"data/${sname}_ts_eb_real\"
\"data/${sname}_coll\" \"data/${sname}_qa\" \"${setup}\" -1
- \"data/${sname}_sign\" \"data/${sname}_beam\" )
+ \"data/${sname}_sign\" \"data/${sname}_beam\" \"qa_config.cbm.yaml\")
set_tests_properties(${testname} PROPERTIES
TIMEOUT ${timeOutTime}
FAIL_REGULAR_EXPRESSION "segmentation violation"
@@ -481,7 +486,7 @@ EndIf() # If(DEFINED ENV{RAW_DATA_PATH} )
#####################
# ============================================================================
-Install(FILES .rootrc run_tra_file.C run_tra_beam.C run_digi.C run_reco.C run_unpack_online.C run_unpack_tsa.C
+Install(FILES .rootrc run_tra_file.C run_tra_beam.C run_digi.C run_reco.C run_unpack_online.C run_unpack_tsa.C qa_config.cbm.yaml
DESTINATION share/cbmroot/macro/run
)
Install(PROGRAMS run_tests.sh
diff --git a/macro/run/qa_config.cbm.yaml b/macro/run/qa_config.cbm.yaml
new file mode 100644
index 0000000000000000000000000000000000000000..6d23fa87225a98f1c8e9d22bd8235b746d144a8d
--- /dev/null
+++ b/macro/run/qa_config.cbm.yaml
@@ -0,0 +1,288 @@
+#
+# -- Configuration for QA tasks --
+#
+# Configuration file defines histogram properties for different QA
+# tasks in Full-CBM setup.
+#
+# Properties of histograms:
+#
+# TH1:
+# title - title of histogram
+# nbinsx - number of bins in the X axis
+# minx - lower bound of the X axis
+# maxx - upper bound of the X axis
+#
+# TH2:
+# title - title of histogram
+# nbinsx - number of bins in the X axis
+# minx - lower bound of the X axis
+# maxx - upper bound of the X axis
+# nbinsy - number of bins in the Y axis
+# miny - lower bound of the Y axis
+# maxy - upper bound of the Y axis
+#
+# TH3:
+# title - title of histogram
+# nbinsx - number of bins in the X axis
+# minx - lower bound of the X axis
+# maxx - upper bound of the X axis
+# nbinsy - number of bins in the Y axis
+# miny - lower bound of the Y axis
+# maxy - upper bound of the Y axis
+# nbinsz - number of bins in the Z axis
+# minz - lower bound of the Z axis
+# maxz - upper bound of the Z axis
+#
+# TProfile:
+# title - title of the profile
+# nbinsx - number of bins in the X axis
+# minx - lower bound of the X axis
+# maxx - upper bound of the X axis
+# miny - lower bound of the Y axis (optional)
+# maxy - upper bound of the Y axis (optional)
+# opt - option (optional)
+#
+# TProfile2D:
+# title - title of the profile
+# nbinsx - number of bins in the X axis
+# minx - lower bound of the X axis
+# maxx - upper bound of the X axis
+# nbinsy - number of bins in the Y axis
+# miny - lower bound of the Y axis
+# maxy - upper bound of the Y axis
+# minz - lower bound of the Z axis (optional)
+# maxz - upper bound of the Z axis (optional)
+# opt - option (optional)
+
+qa:
+ CbmCaOutputQa:
+ histograms:
+ #
+ # Reconstructed tracks distributions
+ #
+ reco_p:
+ title: "Momentum of reconstructed track;(p/q)_{reco} [GeV/ec]"
+ nbinsx: 150
+ minx: 0.
+ maxx: 15.
+ reco_pt:
+ title: "Transverse momentum of reconstructed track;(p_{T}/q)_{reco} [GeV/ec]"
+ nbinsx: 100
+ minx: 0.
+ maxx: 10.
+ reco_phi:
+ title: "Azimuthal angle of reconstructed track;#phi_{reco} [rad]"
+ nbinsx: 68
+ minx: -3.2
+ maxx: +3.2
+ reco_tx:
+ title: "Slope along x-axis of reconstructed tracks;t_{x}"
+ nbinsx: 50
+ minx: -2.
+ maxx: +2.
+ reco_ty:
+ title: "Slope along y-axis of reconstructed tracks;t_{y}"
+ nbinsx: 50
+ minx: -2.
+ maxx: +2.
+ reco_fhitR:
+ title: "Distance of the first hit from z-axis of reconstructed track;R [cm]"
+ nbinsx: 50
+ minx: 0.
+ maxx: 150.
+ reco_nhits:
+ title: "Number of hits in a reconstructed track;N_{hits}"
+ nbinsx: 10
+ minx: -0.5
+ maxx: 9.5
+ reco_fsta:
+ title: "First station of reconstructed track;ID_{f.st.}"
+ nbinsx: 15
+ minx: -0.5
+ maxx: 14.5
+ reco_purity:
+ title: "Purity of reconstructed track;purity [%]"
+ nbinsx: 15
+ minx: -0.5
+ maxx: 14.5
+ reco_chi2_ndf:
+ title: "Fit #chi^{2}/NDF of reconstructed track;(#chi^{2}/NDF)_{reco}"
+ nbinsx: 50
+ minx: -0.5
+ maxx: 99.5
+ reco_prob:
+ title: "Fit #chi^{2}/NDF of the rest tracks;(#chi^{2}/NDF)_{reco}"
+ nbinsx: 101
+ minx: 0.
+ maxx: 1.01
+ rest_chi2_ndf:
+ title: "Fit probability of reconstructed track;Prob_{reco}"
+ nbinsx: 50
+ minx: -0.5
+ maxx: 99.5
+ rest_prob:
+ title: "Fit probability of the rest tracks;Prob_{reco}"
+ nbinsx: 101
+ minx: 0.
+ maxx: 1.01
+ #
+ # Ghost tracks distributions
+ #
+ ghost_p:
+ title: "Momentum of ghost track;(p/q)_{reco} [GeV/ec]"
+ nbinsx: 150
+ minx: 0. # Gev/ec
+ maxx: 15. # Gev/ec
+ ghost_pt:
+ title: "Transverse momentum of ghost track;(p_{T}/q)_{reco} [GeV/ec]"
+ nbinsx: 100
+ minx: 0.
+ maxx: 10.
+ ghost_phi:
+ title: "Azimuthal angle of ghost track;#phi_{reco} [rad]"
+ nbinsx: 68
+ minx: -3.2
+ maxx: +3.2
+ ghost_tx:
+ title: "Slope along x-axis of ghost tracks;t_{x}"
+ nbinsx: 50
+ minx: -2.
+ maxx: +2.
+ ghost_ty:
+ title: "Slope along y-axis of ghost tracks;t_{y}"
+ nbinsx: 50
+ minx: -2.
+ maxx: +2.
+ ghost_fhitR:
+ title: "Distance of the first hit from z-axis of ghost track;R [cm]"
+ nbinsx: 50
+ minx: 0.
+ maxx: 150.
+ ghost_nhits:
+ title: "Number of hits in a ghost track;N_{hits}"
+ nbinsx: 10
+ minx: -0.5
+ maxx: 9.5
+ ghost_fsta:
+ title: "First station of ghost track;ID_{f.st.}"
+ nbinsx: 15
+ minx: -0.5
+ maxx: 14.5
+ ghost_purity:
+ title: "Purity of ghost track;purity [%]"
+ nbinsx: 15
+ minx: -0.5
+ maxx: 14.5
+ ghost_chi2_ndf:
+ title: "Fit #chi^{2}/NDF of ghost track;(#chi^{2}/NDF)_{reco}"
+ nbinsx: 50
+ minx: -0.5
+ maxx: 99.5
+ ghost_prob:
+ title: "Fit #chi^{2}/NDF of the ghost track;(#chi^{2}/NDF)_{reco}"
+ nbinsx: 101
+ minx: 0.
+ maxx: 1.01
+ ghost_nhits_vs_p:
+ title: "Number of hits vs. momentum for ghost track;(p/q)_{reco} [GeV/ec];N_{hits}"
+ nbinsx: 150
+ minx: 0.
+ maxx: 15.
+ nbinsy: 10
+ miny: -0.5
+ maxy: 9.5
+ ghost_fsta_vs_p:
+ title: "First station vs. momentum for ghost track;(p/q)_{reco} [GeV/ec];ID_{f.st.}"
+ nbinsx: 150
+ minx: 0.
+ maxx: 15.
+ nbinsy: 15
+ miny: -0.5
+ maxy: 14.5
+ ghost_lsta_vs_fsta:
+ title: "Last station vs. first station for ghost track;ID_{first st.};ID_{last st.}"
+ nbinsx: 15
+ minx: -0.5
+ maxx: 14.5
+ nbinsy: 15
+ miny: -0.5
+ maxy: 14.5
+ #
+ # -- Track residuals and pulls --
+ #
+ # Residuals and pulls at first station
+ #
+ fst_res_x:
+ titie: "Residual of x in first track point;(x_{reco} - x_{MC})[cm]"
+ nbinsx: 160
+ minx: -4.
+ maxx: +4.
+ fst_res_y:
+ titie: "Residual of y in first track point;(y_{reco} - y_{MC})[cm]"
+ nbinsx: 180
+ minx: -45.
+ maxx: +45.
+ fst_res_tx:
+ titie: "Residual of t_{x} in first track point;(t^{reco}_{x} - t^{MC}_{x})"
+ nbinsx: 80
+ minx: -4.
+ maxx: +4.
+ fst_res_ty:
+ titie: "Residual of t_{y} in first track point;(t^{reco}_{y} - t^{MC}_{y})"
+ nbinsx: 80
+ minx: -4.
+ maxx: +4.
+ fst_res_qp:
+ titie: "Residual of q/p in first track point;((q/p)_{reco} - (q/p)_{MC}) [GeV/ec]"
+ nbinsx: 150
+ minx: -15.
+ maxx: +15.
+ fst_res_time:
+ titie: "Residual of time in first track point;((t)_{reco} - (t)_{MC}) [ns]"
+ nbinsx: 100
+ minx: -10.
+ maxx: +10.
+ fst_res_v:
+ titie: "Residual of velocity in first track point;(v_{reco} - v_{MC}) [c]"
+ nbinsx: 100
+ minx: -10.
+ maxx: +10.
+ fst_pull_x:
+ titie: "Pull of x in first track point;(x_{reco} - x_{MC}) / #sigma_{x}"
+ nbinsx: 120
+ minx: -6.
+ maxx: +6.
+ fst_pull_y:
+ titie: "Pull of y in first track point;(y_{reco} - y_{MC}) / #sigma_{y}"
+ nbinsx: 120
+ minx: -6.
+ maxx: +6.
+ fst_pull_tx:
+ titie: "Pull of t_{x} in first track point;(t^{reco}_{x} - t^{MC}_{x}) / #sigma_{t_{x}}"
+ nbinsx: 120
+ minx: -6.
+ maxx: +6.
+ fst_pull_ty:
+ titie: "Pull of t_{y} in first track point;(t^{reco}_{y} - t^{MC}_{y}) / #sigma_{t_{y}}"
+ nbinsx: 120
+ minx: -6.
+ maxx: +6.
+ fst_pull_qp:
+ titie: "Pull of q/p in first track point;((q/p)_{reco} - (q/p)_{MC}) / #sigma_{q/p}"
+ nbinsx: 120
+ minx: -6.
+ maxx: +6.
+ fst_pull_time:
+ titie: "Pull of time in first track point;((t)_{reco} - (t)_{MC}) / #sigma_{t}"
+ nbinsx: 120
+ minx: -6.
+ maxx: +6.
+ fst_pull_v:
+ titie: "Pull of velocity in first track point;(v_{reco} - v_{MC}) / #sigma_{v}"
+ nbinsx: 120
+ minx: -6.
+ maxx: +6.
+
+
+
+
diff --git a/macro/run/run_qa.C b/macro/run/run_qa.C
index 5a2d18f6edf44b521ded5b1ff8b31b80b476c60f..bb9f10825399cf0bb53351f01f8195af83c09737 100644
--- a/macro/run/run_qa.C
+++ b/macro/run/run_qa.C
@@ -18,17 +18,13 @@
// Includes needed for IDE
#if !defined(__CLING__)
-#include "CbmCaInputQaMuch.h"
-#include "CbmCaInputQaSts.h"
-#include "CbmCaInputQaTof.h"
-#include "CbmCaInputQaTrd.h"
-#include "CbmCaOutputQa.h"
#include "CbmDefs.h"
#include "CbmMCDataManager.h"
#include "CbmMuchDigitizerQa.h"
#include "CbmMuchHitFinderQa.h"
#include "CbmMuchTransportQa.h"
#include "CbmSetup.h"
+#include "CbmStsFindTracksQa.h"
#include <FairFileSource.h>
#include <FairMonitor.h>
@@ -45,7 +41,7 @@
void run_qa(TString dataTra = "data/sis100_muon_jpsi_test", TString dataRaw = "data/sis100_muon_jpsi_test",
TString dataReco = "data/sis100_muon_jpsi_test", TString dataPar = "data/sis100_muon_jpsi_test",
TString dataSink = "data/sis100_muon_jpsi_test", TString setup = "sis100_muon_jpsi", Int_t nEvents = -1,
- TString dataTra2 = "", TString dataTra3 = "")
+ TString dataTra2 = "", TString dataTra3 = "", const char* configName = "./qa_config.fcbm.yaml")
{
gROOT->SetBatch(kTRUE);
@@ -58,9 +54,12 @@ void run_qa(TString dataTra = "data/sis100_muon_jpsi_test", TString dataRaw = "d
fair::Logger::DefineVerbosity(
"user1", fair::VerbositySpec::Make(fair::VerbositySpec::Info::severity, fair::VerbositySpec::Info::file_line));
FairLogger::GetLogger()->SetLogVerbosityLevel("user1");
+ FairLogger::GetLogger()->SetColoredLog(true);
// ------------------------------------------------------------------------
// ----- Environment --------------------------------------------------
+ bool bUseMC = true; // MC flag: used or not
+ int verbose = 3;
TString myName = "run_qa"; // this macro's name for screen output
TString srcDir = gSystem->Getenv("VMCWORKDIR"); // top source directory
// ------------------------------------------------------------------------
@@ -78,8 +77,9 @@ void run_qa(TString dataTra = "data/sis100_muon_jpsi_test", TString dataRaw = "d
// ----- Load the geometry setup -------------------------------------
std::cout << std::endl;
std::cout << "-I- " << myName << ": Loading setup " << setup << std::endl;
- CbmSetup* geo = CbmSetup::Instance();
+ auto* geo = CbmSetup::Instance();
geo->LoadSetup(setup);
+
// You can modify the pre-defined setup by using
// CbmSetup::Instance()->RemoveModule(ESystemId) or
// CbmSetup::Instance()->SetModule(ESystemId, const char*, Bool_t) or
@@ -96,6 +96,13 @@ void run_qa(TString dataTra = "data/sis100_muon_jpsi_test", TString dataRaw = "d
bool bUseTrd = geo->IsActive(ECbmModuleId::kTrd);
bool bUseTof = geo->IsActive(ECbmModuleId::kTof);
bool bUsePsd = geo->IsActive(ECbmModuleId::kPsd);
+ std::cout << " MVD: " << (bUseMvd ? "ON" : "OFF") << '\n';
+ std::cout << " STS: " << (bUseSts ? "ON" : "OFF") << '\n';
+ std::cout << " RICH: " << (bUseRich ? "ON" : "OFF") << '\n';
+ std::cout << " MUCH: " << (bUseMuch ? "ON" : "OFF") << '\n';
+ std::cout << " TRD: " << (bUseTrd ? "ON" : "OFF") << '\n';
+ std::cout << " TOF: " << (bUseTof ? "ON" : "OFF") << '\n';
+ std::cout << " PSD: " << (bUsePsd ? "ON" : "OFF") << '\n';
// ------------------------------------------------------------------------
// ----- Parameter files as input to the runtime database -------------
@@ -173,6 +180,7 @@ void run_qa(TString dataTra = "data/sis100_muon_jpsi_test", TString dataRaw = "d
mcManager->AddFile(traFile);
if (!dataTra2.IsNull()) mcManager->AddFile(tra2File);
if (!dataTra3.IsNull()) mcManager->AddFile(tra3File);
+
run->AddTask(mcManager);
// ------------------------------------------------------------------------
@@ -189,6 +197,10 @@ void run_qa(TString dataTra = "data/sis100_muon_jpsi_test", TString dataRaw = "d
CbmMuchHitFinderQa* muchHitFinderQa = new CbmMuchHitFinderQa();
muchHitFinderQa->SetGeoFileName(muchParFile);
run->AddTask(muchHitFinderQa);
+
+ auto* pCaInputQaMuch = new CbmCaInputQaMuch(verbose, bUseMC);
+ pCaInputQaMuch->SetEfficiencyThrsh(0.5, 0, 100);
+ run->AddTask(pCaInputQaMuch);
}
// ----- TRD QA ---------------------------------
@@ -206,36 +218,43 @@ void run_qa(TString dataTra = "data/sis100_muon_jpsi_test", TString dataRaw = "d
run->AddTask(trdHitProducerQa);
run->AddTask(new CbmTrdCalibTracker());
run->AddTask(new CbmTrackerInputQaTrd()); // Tracker requirements to TRD
+ auto* pCaInputQaTrd = new CbmCaInputQaTrd(verbose, bUseMC);
+ pCaInputQaTrd->SetEfficiencyThrsh(0.5, 0, 100);
+ run->AddTask(pCaInputQaTrd);
}
// ------------------------------------------------------------------------
// ----- TOF QA ---------------------------------
if (CbmSetup::Instance()->IsActive(ECbmModuleId::kTof)) {
//run->AddTask(new CbmTrackerInputQaTof()); // Tracker requirements to TOF
+ auto* pCaInputQaTof = new CbmCaInputQaTof(verbose, bUseMC);
+ pCaInputQaTof->SetEfficiencyThrsh(0.5, 0, 100);
+ run->AddTask(pCaInputQaTof);
}
// ------------------------------------------------------------------------
// ----- STS QA ---------------------------------
if (CbmSetup::Instance()->IsActive(ECbmModuleId::kSts)) {
//run->AddTask(new CbmStsDigitizeQa()); //opens lots of windows
- //run->AddTask(new CbmStsFindTracksQa());
- //run->AddTask(new CbmTrackingInputQaSts());
+ run->AddTask(new CbmStsFindTracksQa());
+ auto* pCaInputQaSts = new CbmCaInputQaSts(verbose, bUseMC);
+ pCaInputQaSts->SetEfficiencyThrsh(0.5, 0, 100);
+ run->AddTask(pCaInputQaSts);
}
// ------------------------------------------------------------------------
- // ----- Event builder QA ------------------------------------------------
+ // ----- Event builder QA ---------------------------------
CbmBuildEventsQa* evBuildQA = new CbmBuildEventsQa();
- //run->AddTask(evBuildQA);
+ run->AddTask(evBuildQA);
// ------------------------------------------------------------------------
// ----- Tracking QA ------------------------------------------------------
- int verbose = 3;
- bool isMCUsed = true;
TString caParFile = recFile;
caParFile.ReplaceAll(".root", ".L1Parameters.dat");
- auto* pCaOutputQa = new cbm::ca::OutputQa(verbose, isMCUsed);
+ auto* pCaOutputQa = new cbm::ca::OutputQa(verbose, bUseMC);
pCaOutputQa->SetStsTrackingMode();
+ pCaOutputQa->SetConfigName(configName);
pCaOutputQa->ReadParameters(caParFile.Data());
pCaOutputQa->SetUseMvd(bUseMvd);
pCaOutputQa->SetUseSts(bUseSts);
diff --git a/reco/L1/CbmCaMCModule.cxx b/reco/L1/CbmCaMCModule.cxx
index 126cc20bda252e963c7e4fe1645bd6cdb2f55663..3f56f34819c291638755e228493413d5d906e5c6 100644
--- a/reco/L1/CbmCaMCModule.cxx
+++ b/reco/L1/CbmCaMCModule.cxx
@@ -115,13 +115,6 @@ try {
LOG(info) << "CA MC Module: initializing CA tracking Monte-Carlo module... \033[1;32mDone!\033[0m";
- LOG(info) << "Detector subsystems used:";
- LOG(info) << "\tMVD: " << (fbUseMvd ? "ON" : "OFF");
- LOG(info) << "\tSTS: " << (fbUseSts ? "ON" : "OFF");
- LOG(info) << "\tMuCh: " << (fbUseMuch ? "ON" : "OFF");
- LOG(info) << "\tTRD: " << (fbUseTrd ? "ON" : "OFF");
- LOG(info) << "\tTOF: " << (fbUseTof ? "ON" : "OFF");
-
return true;
}
catch (const std::logic_error& error) {
@@ -134,9 +127,7 @@ catch (const std::logic_error& error) {
//
void CbmCaMCModule::InitEvent(CbmEvent* /*pEvent*/)
{
- std::cout << "\033[1;32mCALL CbmCAMCModule::InitEvent\033[0m\n";
-
- // ------ Fill a set of file and event indexes
+ // Fill a set of file and event indexes
fFileEventIDs.clear();
if (fbLegacyEventMode) {
int iFile = FairRunAna::Instance()->GetEventHeader()->GetInputFileId();
@@ -152,18 +143,19 @@ void CbmCaMCModule::InitEvent(CbmEvent* /*pEvent*/)
}
}
- // ------ Read data
+ // Read data
fpMCData->Clear();
this->ReadMCTracks();
this->ReadMCPoints();
- L1_SHOW(fpMCData->GetNofTracks());
- L1_SHOW(fpMCData->GetNofPoints());
- // ------ Prepare tracks: set point indexes and redefine indexes from external to internal containers
+ // Prepare tracks: set point indexes and redefine indexes from external to internal containers
for (auto& aTrk : fpMCData->GetTrackContainer()) {
aTrk.SortPointIndexes(
[&](const int& iPl, const int& iPr) { return fpMCData->GetPoint(iPl).GetZ() < fpMCData->GetPoint(iPr).GetZ(); });
}
+
+ // ------ Match reconstructed and MC data
+ this->MatchRecoAndMC();
}
// ---------------------------------------------------------------------------------------------------------------------
@@ -172,16 +164,14 @@ void CbmCaMCModule::ProcessEvent(CbmEvent*) {}
// ---------------------------------------------------------------------------------------------------------------------
//
-void CbmCaMCModule::InitTrackInfo(const L1Vector<CbmL1HitDebugInfo>& vHits)
+void CbmCaMCModule::InitTrackInfo()
{
- LOG(info) << "\033[1;32m!!!! FLAG 1\033[0m";
// ----- Initialize stations arrangement and hit indexes
- fpMCData->InitTrackInfo(vHits);
+ fpMCData->InitTrackInfo(*fpvQaHits);
- LOG(info) << "\033[1;32m!!!! FLAG 2\033[0m";
// ----- Define reconstructable and additional flags
for (auto& aTrk : fpMCData->GetTrackContainer()) {
- bool isRec = true; // is track reconstructable
+ bool isRec = true; // the track can be reconstructed
// Cut on momentum
isRec &= aTrk.GetP() > CbmL1Constants::MinRecoMom;
@@ -213,7 +203,6 @@ void CbmCaMCModule::InitTrackInfo(const L1Vector<CbmL1HitDebugInfo>& vHits)
aTrk.SetFlagReconstructable(isRec);
aTrk.SetFlagAdditional(isAdd);
}
- LOG(info) << "\033[1;32m!!!! FLAG 3\033[0m";
}
// ---------------------------------------------------------------------------------------------------------------------
@@ -227,32 +216,21 @@ void CbmCaMCModule::Finish() {}
// ---------------------------------------------------------------------------------------------------------------------
//
-void CbmCaMCModule::MatchPointsWithHits()
+void CbmCaMCModule::MatchRecoAndMC()
{
- // FIXME: Cleaning up makes it safer, but probably is not needed after old code removal
- for (auto& hit : (*fpvHitIds)) {
- hit.mcPointIds.clear();
- }
-
- int nHits = fpvHitIds->size();
- for (int iH = 0; iH < nHits; ++iH) {
- auto& hit = (*fpvHitIds)[iH];
- int iP = fpMCData->GetPointIndexOfHit(iH);
- if (iP > -1) {
- hit.mcPointIds.push_back_no_warning(iP);
- fpMCData->GetPoint(iP).AddHitID(iH);
- }
- } // loop over hit indexes: end
+ this->MatchPointsAndHits();
+ this->MatchRecoAndMCTracks();
+ this->InitTrackInfo();
}
+
// ---------------------------------------------------------------------------------------------------------------------
//
template<>
-int CbmCaMCModule::MatchHitWithMc<L1DetectorID::kMvd>(int iHit) const
+int CbmCaMCModule::MatchHitWithMc<L1DetectorID::kMvd>(int iHitExt) const
{
int iPoint = -1;
if (fpMvdHitMatches) {
- int iHitExt = -(1 + iHit); // TODO: SZh 28.08.2022: this should be replaced with iHitExt = hit.extIdex
const auto* hitMatch = dynamic_cast<CbmMatch*>(fpMvdHitMatches->At(iHitExt));
assert(hitMatch);
if (hitMatch->GetNofLinks() > 0 && hitMatch->GetLink(0).GetIndex() < fvNofPointsUsed[int(L1DetectorID::kMvd)]) {
@@ -265,10 +243,10 @@ int CbmCaMCModule::MatchHitWithMc<L1DetectorID::kMvd>(int iHit) const
// ---------------------------------------------------------------------------------------------------------------------
//
template<>
-int CbmCaMCModule::MatchHitWithMc<L1DetectorID::kSts>(int iHit) const
+int CbmCaMCModule::MatchHitWithMc<L1DetectorID::kSts>(int iHitExt) const
{
int iPoint = -1;
- const auto* sh = dynamic_cast<CbmStsHit*>(fpStsHits->At(iHit));
+ const auto* sh = dynamic_cast<CbmStsHit*>(fpStsHits->At(iHitExt));
// Match MC point
if (fpStsClusterMatches) {
@@ -321,10 +299,10 @@ int CbmCaMCModule::MatchHitWithMc<L1DetectorID::kSts>(int iHit) const
// ---------------------------------------------------------------------------------------------------------------------
//
template<>
-int CbmCaMCModule::MatchHitWithMc<L1DetectorID::kMuch>(int iHit) const
+int CbmCaMCModule::MatchHitWithMc<L1DetectorID::kMuch>(int iHitExt) const
{
int iPoint = -1;
- const auto* hitMatchMuch = dynamic_cast<CbmMatch*>(fpMuchHitMatches->At(iHit));
+ const auto* hitMatchMuch = dynamic_cast<CbmMatch*>(fpMuchHitMatches->At(iHitExt));
if (hitMatchMuch) {
for (int iLink = 0; iLink < hitMatchMuch->GetNofLinks(); ++iLink) {
if (hitMatchMuch->GetLink(iLink).GetIndex() < fvNofPointsUsed[int(L1DetectorID::kMuch)]) {
@@ -346,10 +324,10 @@ int CbmCaMCModule::MatchHitWithMc<L1DetectorID::kMuch>(int iHit) const
// ---------------------------------------------------------------------------------------------------------------------
//
template<>
-int CbmCaMCModule::MatchHitWithMc<L1DetectorID::kTrd>(int iHit) const
+int CbmCaMCModule::MatchHitWithMc<L1DetectorID::kTrd>(int iHitExt) const
{
int iPoint = -1;
- const auto* hitMatch = dynamic_cast<const CbmMatch*>(fpTrdHitMatches->At(iHit));
+ const auto* hitMatch = dynamic_cast<const CbmMatch*>(fpTrdHitMatches->At(iHitExt));
if (hitMatch) {
int iMC = -1;
if (hitMatch->GetNofLinks() > 0) {
@@ -371,10 +349,10 @@ int CbmCaMCModule::MatchHitWithMc<L1DetectorID::kTrd>(int iHit) const
// ---------------------------------------------------------------------------------------------------------------------
//
template<>
-int CbmCaMCModule::MatchHitWithMc<L1DetectorID::kTof>(int iHit) const
+int CbmCaMCModule::MatchHitWithMc<L1DetectorID::kTof>(int iHitExt) const
{
int iPoint = -1;
- const auto* hitMatch = dynamic_cast<const CbmMatch*>(fpTofHitMatches->At(iHit));
+ const auto* hitMatch = dynamic_cast<const CbmMatch*>(fpTofHitMatches->At(iHitExt));
if (hitMatch) {
for (int iLink = 0; iLink < hitMatch->GetNofLinks(); ++iLink) {
int iMc = hitMatch->GetLink(iLink).GetIndex();
@@ -391,6 +369,55 @@ int CbmCaMCModule::MatchHitWithMc<L1DetectorID::kTof>(int iHit) const
return iPoint;
}
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void CbmCaMCModule::MatchPointsAndHits()
+{
+ // Clear point indexes of each hit
+ // NOTE: at the moment only the best point is stored to the hit, but we keep the possibility for keeping many
+ // points for possible future module extension.
+ for (auto& hit : (*fpvQaHits)) {
+ hit.mcPointIds.clear();
+ }
+
+ // Match hit with point for different detectors
+ auto DoMatch = [&](int iH, int iP) constexpr
+ {
+ if (iP > -1) {
+ auto& hit = (*fpvQaHits)[iH];
+ hit.mcPointIds.push_back_no_warning(iP);
+ fpMCData->GetPoint(iP).AddHitID(iH);
+ }
+ };
+
+ int iH = 0;
+ if (fbUseMvd) {
+ for (iH = (*fpvFstHitId)[int(L1DetectorID::kMvd)]; iH < (*fpvFstHitId)[int(L1DetectorID::kMvd) + 1]; ++iH) {
+ DoMatch(iH, MatchHitWithMc<L1DetectorID::kMvd>((*fpvQaHits)[iH].ExtIndex));
+ }
+ }
+ if (fbUseSts) {
+ for (iH = (*fpvFstHitId)[int(L1DetectorID::kSts)]; iH < (*fpvFstHitId)[int(L1DetectorID::kSts) + 1]; ++iH) {
+ DoMatch(iH, MatchHitWithMc<L1DetectorID::kSts>((*fpvQaHits)[iH].ExtIndex));
+ }
+ }
+ if (fbUseSts) {
+ for (iH = (*fpvFstHitId)[int(L1DetectorID::kMuch)]; iH < (*fpvFstHitId)[int(L1DetectorID::kMuch) + 1]; ++iH) {
+ DoMatch(iH, MatchHitWithMc<L1DetectorID::kMuch>((*fpvQaHits)[iH].ExtIndex));
+ }
+ }
+ if (fbUseTrd) {
+ for (iH = (*fpvFstHitId)[int(L1DetectorID::kTrd)]; iH < (*fpvFstHitId)[int(L1DetectorID::kTrd) + 1]; ++iH) {
+ DoMatch(iH, MatchHitWithMc<L1DetectorID::kTrd>((*fpvQaHits)[iH].ExtIndex));
+ }
+ }
+ if (fbUseTof) {
+ for (iH = (*fpvFstHitId)[int(L1DetectorID::kTof)]; iH < (*fpvFstHitId)[int(L1DetectorID::kTof) + 1]; ++iH) {
+ DoMatch(iH, MatchHitWithMc<L1DetectorID::kTof>((*fpvQaHits)[iH].ExtIndex));
+ }
+ }
+}
+
// ---------------------------------------------------------------------------------------------------------------------
//
void CbmCaMCModule::MatchRecoAndMCTracks()
@@ -402,7 +429,7 @@ void CbmCaMCModule::MatchRecoAndMCTracks()
auto& mNofHitsVsMCTrkID = trkRe.hitMap;
mNofHitsVsMCTrkID.clear();
for (int iH : trkRe.Hits) {
- for (int iP : (*fpvHitIds)[iH].mcPointIds) {
+ for (int iP : (*fpvQaHits)[iH].mcPointIds) {
assert(iP > -1); // Should never be triggered
int iTmc = fpMCData->GetPoint(iP).GetTrackId();
if (mNofHitsVsMCTrkID.find(iTmc) == mNofHitsVsMCTrkID.end()) { mNofHitsVsMCTrkID[iTmc] = 1; }
@@ -472,20 +499,18 @@ void CbmCaMCModule::SetDetector(L1DetectorID detID, bool flag)
//
void CbmCaMCModule::CheckInit() const
{
- // ----- Check parameters
+ // Check parameters
if (!fpParameters.get()) { throw std::logic_error("Tracking parameters object was not defined"); }
+ // Check output data containers
if (!fpMCData) { throw std::logic_error("MC data object was not registered"); }
-
if (!fpvRecoTracks) { throw std::logic_error("Reconstructed track container was not registered"); }
-
if (!fpvHitIds) { throw std::logic_error("Hit index container was not registered"); }
-
- if (!fpInputData) { throw std::logic_error("Input data object was not registered"); }
+ if (!fpvQaHits) { throw std::logic_error("QA hit container was not registered"); }
+ if (!fpvFstHitId) { throw std::logic_error("Array of first hit indexes in each detector was not registered"); }
// Check event list
if (!fbLegacyEventMode && !fpMCEventList) { throw std::logic_error("MC event list was not found"); }
-
if (!fbLegacyEventMode && !fpTimeSlice) { throw std::logic_error("Time slice was not found"); }
// Tracks branch
@@ -716,7 +741,7 @@ void CbmCaMCModule::ReadMCTracks()
else {
// This is a secondary track, mother ID should be recalculated for the internal array
int motherId = fpMCData->FindInternalTrackIndex(extMotherId, iEvent, iFile);
- assert(motherId > -1);
+ if (motherId == -1) { motherId = -3; } // Mother is neutral particle, which is rejected
aTrk.SetMotherId(motherId);
}
diff --git a/reco/L1/CbmCaMCModule.h b/reco/L1/CbmCaMCModule.h
index 224fa600314c830755a22794ddb68e9bffc633a2..61efa6695f364d21096c262bbe5e5d9428cbf075 100644
--- a/reco/L1/CbmCaMCModule.h
+++ b/reco/L1/CbmCaMCModule.h
@@ -11,6 +11,7 @@
#define CbmCaMCModule_h 1
#include "CbmL1DetectorID.h"
+#include "CbmL1Hit.h"
#include "CbmL1Track.h"
#include "CbmLink.h"
#include "CbmMCDataArray.h"
@@ -95,41 +96,28 @@ public:
bool InitRun();
/// @brief Initializes MC track
- /// @param vHits Vector of hit objects
///
/// Initializes information about arrangement of points and hits of MC tracks within stations and the status
- /// of track reconstructability, calculates max number of points and hits on a station, number of consecutive
- /// stations containing a hit or point and number of stations and points with hits. Provides the determination
- /// of track reconstructability status.
- void InitTrackInfo(const L1Vector<CbmL1HitDebugInfo>& vHits);
+ /// of track ability to be reconstructed, calculates max number of points and hits on a station, number of
+ /// consecutive stations containing a hit or point and number of stations and points with hits.
+ void InitTrackInfo();
+
+ /// @brief Match reconstructed and MC data
+ ///
+ /// Runs matching of hits with points and reconstructed tracks with
+ void MatchRecoAndMC();
/// @brief Processes event
///
/// Fills histograms and tables, should be called after the tracking done
void ProcessEvent(CbmEvent* pEvent);
- /// @brief Matches hit with MC point
- /// @tparam DetId Detector ID
- /// @param iHit External index of hit
- /// @return MC-point index in fvMCPoints array
- ///
- /// This method finds a match for a given hit or matches for hits clusters (in case of STS), finds the best
- /// link in the match and returns the corresponding global ID of the MC points
- template<L1DetectorID DetId>
- int MatchHitWithMc(int iHit) const;
-
- /// @brief Assigns MC point indexes to each hit and hit indexes to each MC point
- void MatchPointsWithHits();
-
- /// @brief Matches reconstructed tracks with MC tracks
- ///
- /// In the procedure, the maps of associated MC track indexes to corresponding number of hits are filled out and the
- /// max purity is calculated for each reconstructed track in the TS/event. If the value of purity for a given MC track
- /// is larger then a threshold, the corresponding track index is saved to the reconstructed track object, in the same
- /// time the index of the reconstructed track object is saved to this MC track object. If purity for the MC track does
- /// not pass the threshold, its index will not be accounted as a matched track, and the index of reconstructed track
- /// will be added as an index of "touched" track.
- void MatchRecoAndMCTracks();
+ /// @brief Sets first hit indexes container in different detectors
+ /// @param source Array of indexes
+ void RegisterFirstHitIndexes(const std::array<int, L1Constants::size::kMaxNdetectors + 1>& source)
+ {
+ fpvFstHitId = &source;
+ }
/// @brief Sets used detector subsystems
/// @param detID Id of detector
@@ -153,28 +141,23 @@ public:
/// @brief Registers hit index container
/// @param vHitIds Reference to hit index container
- void RegisterHitIndexContainer(L1Vector<CbmL1HitDebugInfo>& vHitIds) { fpvHitIds = &vHitIds; }
-
- /// @brief Registers input data object
- /// @param inputData Instance of the input data object
- void RegisterInputData(L1InputData& inputData) { fpInputData = &inputData; }
-
- // /// @brief Registers debug hit container
- // /// @param vDbgHits Reference to debug hit container
- // void RegisterHitDebugInfoContainer(L1Vector<CbmL1HitStore>& vDbgHits) { fpvDbgHits = &vDbgHits; }
+ void RegisterHitIndexContainer(L1Vector<CbmL1HitId>& vHitIds) { fpvHitIds = &vHitIds; }
/// @brief Registers CA parameters object
/// @param pParameters A shared pointer to the parameters object
void RegisterParameters(std::shared_ptr<L1Parameters>& pParameters) { fpParameters = pParameters; }
+ /// @brief Registers debug hit container
+ /// @param vQaHits Reference to debug hit container
+ void RegisterQaHitContainer(L1Vector<CbmL1HitDebugInfo>& vQaHits) { fpvQaHits = &vQaHits; }
private:
/// @brief Calculates global index of MC point
/// @param iPointLocal Local index of MC point
/// @param detID Detector ID
///
- /// Calculates global index of MC point as a sum of a given local index and total provided number of points in previous
- /// detector subsystem
+ /// The function calculates global index of MC point as a sum of a given local index and total provided number of
+ /// points in previous detector subsystem.
int CalcGlobMCPointIndex(int iPointLocal, L1DetectorID detID) const
{
return iPointLocal + std::accumulate(fvNofPointsOrig.cbegin(), fvNofPointsOrig.cbegin() + int(detID), 0);
@@ -184,6 +167,31 @@ private:
/// @note The function throws std::logic_error, if initialization is incomplete
void CheckInit() const;
+ /// @brief Matches hit with MC point
+ /// @tparam DetId Detector ID
+ /// @param iHitExt External index of hit
+ /// @return MC-point index in fvMCPoints array
+ ///
+ /// This method finds a match for a given hit or matches for hits clusters (in case of STS), finds the best
+ /// link in the match and returns the corresponding global ID of the MC points.
+ template<L1DetectorID DetId>
+ int MatchHitWithMc(int iHitExt) const;
+
+ /// @brief Match MC points and reconstructed hits
+ ///
+ /// Writes indexes of MC points to each hit and indexes of hits to each MC point.
+ void MatchPointsAndHits();
+
+ /// @brief Matches reconstructed tracks with MC tracks
+ ///
+ /// In the procedure, the maps of associated MC track indexes to corresponding number of hits are filled out and the
+ /// max purity is calculated for each reconstructed track in the TS/event. If the value of purity for a given MC track
+ /// is larger then a threshold, the corresponding track index is saved to the reconstructed track object, in the same
+ /// time the index of the reconstructed track object is saved to this MC track object. If purity for the MC track does
+ /// not pass the threshold, its index will not be accounted as a matched track, and the index of reconstructed track
+ /// will be added as an index of "touched" track.
+ void MatchRecoAndMCTracks();
+
/// @brief Reads MC tracks from external trees and saves them to MCDataObject
void ReadMCTracks();
@@ -205,7 +213,6 @@ private:
template<L1DetectorID DetID>
bool FillMCPoint(int iExtId, int iEvent, int iFile, ca::tools::MCPoint& point);
-
std::shared_ptr<L1Parameters> fpParameters = nullptr; ///< Pointer to tracking parameters object
// ------ Flags
@@ -234,8 +241,8 @@ private:
CbmMCDataArray* fpTrdPoints = nullptr; ///< TRD MC-points input container
CbmMCDataArray* fpTofPoints = nullptr; ///< TOF MC-points input container
- std::array<int, L1Constants::size::kMaxNdetectors> fvNofPointsOrig = {}; ///< Number of points by detector provided
- std::array<int, L1Constants::size::kMaxNdetectors> fvNofPointsUsed = {}; ///< Number of points used in performance
+ std::array<int, L1Constants::size::kMaxNdetectors> fvNofPointsOrig = {0}; ///< Number of points by detector provided
+ std::array<int, L1Constants::size::kMaxNdetectors> fvNofPointsUsed = {0}; ///< Number of points used in performance
// Matches
TClonesArray* fpMvdHitMatches = nullptr; ///< Array of MVD hit matches
@@ -250,13 +257,18 @@ private:
// Matching information
std::set<std::pair<int, int>> fFileEventIDs; ///< Set of file and event indexes: first - iFile, second - iEvent
- ca::tools::MCData* fpMCData = nullptr; ///< MC-information in CA tracking internal format (tracks and points)
+ // MC data
+ ca::tools::MCData* fpMCData = nullptr; ///< MC information (hits and tracks) instance
// Reconstructed data
- L1Vector<CbmL1Track>* fpvRecoTracks = nullptr; ///< Pointer to reconstructed track container
- L1Vector<CbmL1HitDebugInfo>* fpvHitIds = nullptr; ///< Pointer to hit index container
- //L1Vector<CbmL1HitStore>* fpvDbgHits = nullptr; ///< Pointer to debug hit container
- L1InputData* fpInputData = nullptr; ///< Pointer to input data object
+ L1Vector<CbmL1Track>* fpvRecoTracks = nullptr; ///< Pointer to reconstructed track container
+ L1Vector<CbmL1HitId>* fpvHitIds = nullptr; ///< Pointer to hit index container
+ L1Vector<CbmL1HitDebugInfo>* fpvQaHits = nullptr; ///< Pointer to QA hit container
+
+ /// @brief Pointer to array of first hit indexes in the detector subsystem
+ ///
+ /// This array must be initialized in the run initialization function.
+ const std::array<int, L1Constants::size::kMaxNdetectors + 1>* fpvFstHitId = nullptr;
};
// **********************************************
diff --git a/reco/L1/CbmCaTimeSliceReader.cxx b/reco/L1/CbmCaTimeSliceReader.cxx
index 8c85c555b1103f71825db762bb45b2dad881efaa..1f55a04df68bee7ba2af52396b25e8661cedcfde 100644
--- a/reco/L1/CbmCaTimeSliceReader.cxx
+++ b/reco/L1/CbmCaTimeSliceReader.cxx
@@ -23,6 +23,7 @@
using cbm::ca::TimeSliceReader;
using L1Constants::clrs::kCL; // clear log
using L1Constants::clrs::kGNb; // green bold log
+using L1Constants::clrs::kRDb; // red bold log
// ---------------------------------------------------------------------------------------------------------------------
//
@@ -41,6 +42,7 @@ void TimeSliceReader::Clear()
// Input branches
fpBrTimeSlice = nullptr;
+ fpParameters = nullptr;
fpBrMvdHits = nullptr;
fpBrStsHits = nullptr;
@@ -56,15 +58,51 @@ void TimeSliceReader::Clear()
// Pointers to output containers
fpvHitIds = nullptr;
- fpvDbgHits = nullptr;
+ fpvQaHits = nullptr;
fpIODataManager = nullptr;
fpvTracks = nullptr;
+
+ fNofHits = 0;
+ fNofHitKeys = 0;
+ fFirstHitKey = 0;
+
+ std::fill(fvHitFirstIndexDet.begin(), fvHitFirstIndexDet.end(), 0);
}
// ---------------------------------------------------------------------------------------------------------------------
//
-void TimeSliceReader::InitRun()
+void TimeSliceReader::CheckInit() const
{
+ // Check parameters
+ if (!fpParameters.get()) { throw std::logic_error("Tracking parameters object was not defined"); }
+ if (!fpvHitIds) { throw std::logic_error("Hit index container was not defined"); }
+
+ if (!fpBrTimeSlice) { throw std::logic_error("Time slice branch is unavailable"); }
+
+ if (fbUseMvd && !fpBrMvdHits) { throw std::logic_error("MVD hits branch is unavailable"); }
+ if (fbUseSts && !fpBrStsHits) { throw std::logic_error("STS hits branch is unavailable"); }
+ if (fbUseMuch && !fpBrMuchHits) { throw std::logic_error("MuCh hits branch is unavailable"); }
+ if (fbUseTrd && !fpBrTrdHits) { throw std::logic_error("TRD hits branch is unavailable"); }
+ if (fbUseTof && !fpBrTofHits) { throw std::logic_error("TOF hits branch is unavailable"); }
+
+ if (fpvTracks) {
+ if (ECbmTrackingMode::kSTS == fTrackingMode) {
+ if (!fpBrRecoTracks) { throw std::logic_error("StsTrack branch is unavailable"); }
+ }
+ else if (ECbmTrackingMode::kMCBM == fTrackingMode) {
+ if (!fpBrRecoTracks) { throw std::logic_error("GlobalTrack branch is unavailable"); }
+ if (fbUseSts && !fpBrStsTracks) { throw std::logic_error("StsTrack branch is unavailable"); }
+ if (fbUseMuch && !fpBrMuchTracks) { throw std::logic_error("MuchTrack branch is unavailable"); }
+ if (fbUseTrd && !fpBrRecoTracks) { throw std::logic_error("TrdTrack branch is unavailable"); }
+ if (fbUseTof && !fpBrRecoTracks) { throw std::logic_error("TofTrack branch is unavailable"); }
+ }
+ }
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+bool TimeSliceReader::InitRun()
+try {
LOG(info) << "TimeSliceReader: initializing run ... ";
// Init tracking detector interfaces
@@ -80,70 +118,48 @@ void TimeSliceReader::InitRun()
LOG_IF(fatal, !pFairManager) << "TimeSliceReader: FairRootManager was not defined";
fpBrTimeSlice = dynamic_cast<CbmTimeSlice*>(pFairManager->GetObject("TimeSlice."));
- LOG_IF(fatal, !fpBrTimeSlice) << "TimeSliceReader: time slice was not defined";
// Init hit branches
- if (fbUseMvd) {
- fpBrMvdHits = dynamic_cast<TClonesArray*>(pFairManager->GetObject("MvdHit"));
- LOG_IF(fatal, !fpBrMvdHits) << "TimeSliceReader: MVD hits were not found";
- }
-
- if (fbUseSts) {
- fpBrStsHits = dynamic_cast<TClonesArray*>(pFairManager->GetObject("StsHit"));
- LOG_IF(fatal, !fpBrStsHits) << "TimeSliceReader: STS hits were not found";
- }
-
- if (fbUseMuch) {
- fpBrMuchHits = dynamic_cast<TClonesArray*>(pFairManager->GetObject("MuchPixelHit"));
- LOG_IF(fatal, !fpBrMuchHits) << "TimeSliceReader: MuCh hits were not found";
- }
-
- if (fbUseTrd) {
- fpBrTrdHits = dynamic_cast<TClonesArray*>(pFairManager->GetObject("TrdHit"));
- LOG_IF(fatal, !fpBrTrdHits) << "TimeSliceReader: TRD hits were not found";
- }
-
- if (fbUseTof) {
- fpBrTofHits = dynamic_cast<TClonesArray*>(pFairManager->GetObject("TofHit"));
- LOG_IF(fatal, !fpBrTofHits) << "TimeSliceReader: TOF hits were not found";
- }
+ if (fbUseMvd) { fpBrMvdHits = dynamic_cast<TClonesArray*>(pFairManager->GetObject("MvdHit")); }
+ if (fbUseSts) { fpBrStsHits = dynamic_cast<TClonesArray*>(pFairManager->GetObject("StsHit")); }
+ if (fbUseMuch) { fpBrMuchHits = dynamic_cast<TClonesArray*>(pFairManager->GetObject("MuchPixelHit")); }
+ if (fbUseTrd) { fpBrTrdHits = dynamic_cast<TClonesArray*>(pFairManager->GetObject("TrdHit")); }
+ if (fbUseTof) { fpBrTofHits = dynamic_cast<TClonesArray*>(pFairManager->GetObject("TofHit")); }
// Init track branches
if (fpvTracks) {
switch (fTrackingMode) {
case ECbmTrackingMode::kSTS:
fpBrRecoTracks = dynamic_cast<TClonesArray*>(pFairManager->GetObject("StsTrack"));
- LOG_IF(fatal, !fpBrRecoTracks) << "TimeSliceReader: StsTrack not found";
break;
case ECbmTrackingMode::kMCBM:
fpBrRecoTracks = dynamic_cast<TClonesArray*>(pFairManager->GetObject("GlobalTrack"));
- LOG_IF(fatal, !fpBrRecoTracks) << "TimeSliceReader: GlobalTrack not found";
- if (fbUseSts) {
- fpBrStsTracks = dynamic_cast<TClonesArray*>(pFairManager->GetObject("StsTrack"));
- LOG_IF(fatal, !fpBrStsTracks) << "TimeSliceReader: StsTrack not found";
- }
- if (fbUseMuch) {
- fpBrMuchTracks = dynamic_cast<TClonesArray*>(pFairManager->GetObject("MuchTrack"));
- LOG_IF(fatal, !fpBrMuchTracks) << "TimeSliceReader: MuchTrack not found";
- }
- if (fbUseTrd) {
- fpBrTrdTracks = dynamic_cast<TClonesArray*>(pFairManager->GetObject("TrdTrack"));
- LOG_IF(fatal, !fpBrTrdTracks) << "TimeSliceReader: TrdTrack not found";
- }
- if (fbUseTof) {
- fpBrTofTracks = dynamic_cast<TClonesArray*>(pFairManager->GetObject("TofTrack"));
- LOG_IF(fatal, !fpBrTofTracks) << "TimeSliceReader: TofTrack not found";
- }
+ if (fbUseSts) { fpBrStsTracks = dynamic_cast<TClonesArray*>(pFairManager->GetObject("StsTrack")); }
+ if (fbUseMuch) { fpBrMuchTracks = dynamic_cast<TClonesArray*>(pFairManager->GetObject("MuchTrack")); }
+ if (fbUseTrd) { fpBrTrdTracks = dynamic_cast<TClonesArray*>(pFairManager->GetObject("TrdTrack")); }
+ if (fbUseTof) { fpBrTofTracks = dynamic_cast<TClonesArray*>(pFairManager->GetObject("TofTrack")); }
break;
}
}
+ // Check initialization
+ this->CheckInit();
+
LOG(info) << "TimeSliceReader: initializing run ... " << kGNb << "done" << kCL;
+ return true;
+}
+catch (const std::logic_error& error) {
+ LOG(info) << "TimeSliceReader: initializing run ... " << kRDb << "failed" << kCL << "\nReason: " << error.what();
+ return false;
}
// ---------------------------------------------------------------------------------------------------------------------
//
-void TimeSliceReader::InitTimeSlice() { LOG(info) << "TimeSliceReader: initializing time slice"; }
+void TimeSliceReader::InitTimeSlice()
+{
+ this->ReadHits();
+ this->ReadRecoTracks();
+}
// ---------------------------------------------------------------------------------------------------------------------
//
@@ -162,21 +178,58 @@ void TimeSliceReader::ReadHits()
if (fbUseTof) { nHitsTot += fpBrTofHits->GetEntriesFast(); }
// Resize the containers
- fpvHitIds->clear();
- fpvHitIds->reserve(nHitsTot);
- fpvDbgHits->clear();
- fpvDbgHits->reserve(nHitsTot);
- fpIODataManager->ResetInputData();
- fpIODataManager->ReserveNhits(nHitsTot);
+ 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) { fNofHits += ReadHitsForDetector<L1DetectorID::kMvd>(fpBrMvdHits); }
- if (fbUseSts) { fNofHits += ReadHitsForDetector<L1DetectorID::kSts>(fpBrStsHits); }
- if (fbUseMuch) { fNofHits += ReadHitsForDetector<L1DetectorID::kMuch>(fpBrMuchHits); }
- if (fbUseTrd) { fNofHits += ReadHitsForDetector<L1DetectorID::kTrd>(fpBrTrdHits); }
- if (fbUseTof) { fNofHits += ReadHitsForDetector<L1DetectorID::kTof>(fpBrTofHits); }
+ 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); }
- fpIODataManager->SetNhitKeys(fNofHitKeys);
+ // Sort debug hits
+ if (fpvQaHits) { this->SortQaHits(); }
}
// ---------------------------------------------------------------------------------------------------------------------
@@ -210,10 +263,12 @@ void TimeSliceReader::ReadRecoTracks()
track.Hits.clear();
track.Hits.reserve(pInputTrack->GetNofHits());
for (int iH = 0; iH < pInputTrack->GetNofMvdHits(); ++iH) {
- track.Hits.push_back(-pInputTrack->GetMvdHitIndex(iH) - 1); // !!!!
- } // iH
+ int iHint = vHitMvdIds[pInputTrack->GetMvdHitIndex(iH)];
+ track.Hits.push_back(iHint); // !!!!
+ } // iH
for (int iH = 0; iH < pInputTrack->GetNofStsHits(); ++iH) {
- track.Hits.push_back(pInputTrack->GetStsHitIndex(iH));
+ int iHint = vHitStsIds[pInputTrack->GetStsHitIndex(iH)];
+ track.Hits.push_back(iHint);
} // iH
} // iT
break;
@@ -225,14 +280,6 @@ void TimeSliceReader::ReadRecoTracks()
}
}
-// ---------------------------------------------------------------------------------------------------------------------
-//
-void TimeSliceReader::RegisterHitDebugInfoContainer(L1Vector<CbmL1HitDebugInfo>& vDbgHits) { fpvDbgHits = &vDbgHits; }
-
-// ---------------------------------------------------------------------------------------------------------------------
-//
-void TimeSliceReader::RegisterHitIndexContainer(L1Vector<CbmL1HitId>& vHitIds) { fpvHitIds = &vHitIds; }
-
// ---------------------------------------------------------------------------------------------------------------------
//
void TimeSliceReader::RegisterIODataManager(std::shared_ptr<L1IODataManager>& pIODataManager)
@@ -241,10 +288,6 @@ void TimeSliceReader::RegisterIODataManager(std::shared_ptr<L1IODataManager>& pI
fpIODataManager = pIODataManager;
}
-// ---------------------------------------------------------------------------------------------------------------------
-//
-void TimeSliceReader::RegisterTracksContainer(L1Vector<CbmL1Track>& vTracks) { fpvTracks = &vTracks; }
-
// ---------------------------------------------------------------------------------------------------------------------
//
void TimeSliceReader::SetDetector(L1DetectorID detID, bool flag)
@@ -260,17 +303,42 @@ void TimeSliceReader::SetDetector(L1DetectorID detID, bool flag)
// ---------------------------------------------------------------------------------------------------------------------
//
-void TimeSliceReader::StoreHitRecord(const HitRecord& hitRecord)
+void TimeSliceReader::SortQaHits()
{
- int iHitGlob = static_cast<int>(fpIODataManager->GetNofHits()); // Current hit global index
+ int nStationsActive = fpParameters->GetNstationsActive();
+ L1Vector<CbmL1HitDebugInfo> vNewHits(fpvQaHits->size());
+ std::vector<int> vHitFirstIndexes(nStationsActive + 1, 0);
+ std::vector<int> vNofHitsStored(nStationsActive, 0);
+
+ // Count number of hits in each station (NOTE: we could use here boarders from the IO data manager, but we would keep
+ // these two procedures independent)
+ for (const auto& hit : (*fpvQaHits)) {
+ ++vHitFirstIndexes[hit.GetStationId() + 1];
+ }
+
+ for (int iSt = 0; iSt < nStationsActive; ++iSt) {
+ vHitFirstIndexes[iSt + 1] += vHitFirstIndexes[iSt];
+ }
+
+ for (const auto& hit : (*fpvQaHits)) {
+ int iSt = hit.GetStationId();
+ vNewHits[vHitFirstIndexes[iSt] + (vNofHitsStored[iSt]++)] = hit;
+ }
+
+ std::swap(vNewHits, (*fpvQaHits));
+}
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void TimeSliceReader::StoreHitRecord(const HitRecord& hitRecord)
+{
// Save the algo hit
if (fpIODataManager.get()) {
L1Hit aHit;
aHit.iSt = hitRecord.fStaId;
aHit.f = hitRecord.fStripF;
aHit.b = hitRecord.fStripB;
- aHit.ID = -1;
+ aHit.ID = static_cast<int>(fpIODataManager->GetNofHits());
aHit.z = hitRecord.fZ;
aHit.u = hitRecord.fU;
aHit.v = hitRecord.fV;
@@ -285,7 +353,7 @@ void TimeSliceReader::StoreHitRecord(const HitRecord& hitRecord)
if (fpvHitIds) { fpvHitIds->emplace_back(hitRecord.fDet, hitRecord.fExtId); }
// Save debug information
- if (fpvDbgHits) {
+ if (fpvQaHits) {
CbmL1HitDebugInfo hitDbg;
hitDbg.Det = hitRecord.fDet;
hitDbg.ExtIndex = hitRecord.fExtId;
@@ -296,6 +364,6 @@ void TimeSliceReader::StoreHitRecord(const HitRecord& hitRecord)
hitDbg.dy = hitRecord.fDy;
hitDbg.dxy = hitRecord.fDxy;
hitDbg.time = hitRecord.fT;
- fpvDbgHits->push_back(hitDbg);
+ fpvQaHits->push_back(hitDbg);
}
}
diff --git a/reco/L1/CbmCaTimeSliceReader.h b/reco/L1/CbmCaTimeSliceReader.h
index 83a4e548e647b25e1fd4c48c67c78d811c8500a6..a92b7809a870f916473f940b56c4fe4a1f78bf23 100644
--- a/reco/L1/CbmCaTimeSliceReader.h
+++ b/reco/L1/CbmCaTimeSliceReader.h
@@ -37,34 +37,35 @@ class L1Parameters;
namespace cbm::ca
{
- /// @brief A structure to store hits information from different detectors in a uniform manner
- struct HitRecord {
- double fX = 0.; ///< x component of hit position [cm]
- double fY = 0.; ///< y component of hit position [cm]
- double fDx = 0.; ///< error of x component of hit position [cm]
- double fDy = 0.; ///< error of y component of hit position [cm]
- double fDxy = 0.; ///< correlation between x and y components [cm]
- double fU = 0.; ///< hit position in direction of front strips [cm]
- double fV = 0.; ///< hit position in direction of back strips [cm]
- double fDu = 0.; ///< hit position error in direction of front strips [cm]
- double fDv = 0.; ///< hit position error in direction of back strips [cm]
- double fZ = 0.; ///< z component of hit position [cm]
- double fT = 0.; ///< time of hit [ns]
- double fDt = 0.; ///< time error of hit [ns]
- int64_t fDataStream = -1; ///< Global index of detector module
- int fExtId = -1; ///< external index of hit
- int fStaId = -1; ///< index of active tracking station
- int fStripF = -1; ///< index of front strip
- int fStripB = -1; ///< index of back strip
- int fDet = -1; ///< detector ID
- };
-
/// @brief A reader of time slice for CA tracker
///
/// The class reads reconstructed hits and reconstructed tracks (optionally) and fills the CA tracking internal
/// data structures.
///
class TimeSliceReader {
+ /// @brief A structure to store hits information from different detectors in a uniform manner
+ ///
+ struct HitRecord {
+ double fX = 0.; ///< x component of hit position [cm]
+ double fY = 0.; ///< y component of hit position [cm]
+ double fDx = 0.; ///< error of x component of hit position [cm]
+ double fDy = 0.; ///< error of y component of hit position [cm]
+ double fDxy = 0.; ///< correlation between x and y components [cm]
+ double fU = 0.; ///< hit position in direction of front strips [cm]
+ double fV = 0.; ///< hit position in direction of back strips [cm]
+ double fDu = 0.; ///< hit position error in direction of front strips [cm]
+ double fDv = 0.; ///< hit position error in direction of back strips [cm]
+ double fZ = 0.; ///< z component of hit position [cm]
+ double fT = 0.; ///< time of hit [ns]
+ double fDt = 0.; ///< time error of hit [ns]
+ int64_t fDataStream = -1; ///< Global index of detector module
+ int fExtId = -1; ///< external index of hit
+ int fStaId = -1; ///< index of active tracking station
+ int fStripF = -1; ///< index of front strip
+ int fStripB = -1; ///< index of back strip
+ int fDet = -1; ///< detector ID
+ };
+
public:
/// @brief Constructor from parameters
/// @param mode Tracking mode
@@ -88,14 +89,31 @@ namespace cbm::ca
/// @brief Clears class content
void Clear();
+ /// @brief Check class initialization
+ /// @note The function throws std::logic_error, if initialization is incomplete
+ void CheckInit() const;
+
+ /// @brief Gets reference to container of first hit indexes in a detector subsystem
+ /// @return Ref. to the container
+ const auto& GetHitFirstIndexDet() const { return fvHitFirstIndexDet; }
+
+ /// @brief Gets number of hits stored for a given detector
+ /// @param iDet Detector ID
+ /// @return Number of hits
+ int GetNofHits(L1DetectorID iDet) const
+ {
+ return fvHitFirstIndexDet[int(iDet) + 1] - fvHitFirstIndexDet[int(iDet)];
+ }
+
/// @brief Run initializer function
+ /// @return Success flag
///
/// Initializes data branches and provides necessary checks in the beginning of the run
- void InitRun();
+ bool InitRun();
/// @brief Reads time slice
///
- /// Reads hits from time slice
+ /// Reads hits and tracks (optionally) from time slice
void InitTimeSlice();
/// @brief Reads hits
@@ -105,20 +123,28 @@ namespace cbm::ca
void ReadRecoTracks();
/// @brief Registers hit debug info container
- /// @param vDbgHits Reference to debug hit container
- void RegisterHitDebugInfoContainer(L1Vector<CbmL1HitDebugInfo>& vDbgHits);
+ /// @param vQaHits Reference to Qa hit container
+ /// @note If no container is registered, all related routines are omitted
+ void RegisterQaHitContainer(L1Vector<CbmL1HitDebugInfo>& vQaHits) { fpvQaHits = &vQaHits; }
/// @brief Registers hit index container
/// @param vHitIds Reference to hits indexes container
- void RegisterHitIndexContainer(L1Vector<CbmL1HitId>& vHitIds);
+ /// @note If no container is registered, all related routines are omitted
+ void RegisterHitIndexContainer(L1Vector<CbmL1HitId>& vHitIds) { fpvHitIds = &vHitIds; }
+
+ /// @brief Registers CA parameters object
+ /// @param pParameters A shared pointer to the parameters object
+ void RegisterParameters(std::shared_ptr<L1Parameters>& pParameters) { fpParameters = pParameters; }
/// @brief Registers the CA IO data manager instance
/// @param pIODataManager Shared pointer to the IO data manager instance
+ /// @note If no container is registered, all related routines are omitted
void RegisterIODataManager(std::shared_ptr<L1IODataManager>& ioDataManager);
/// @brief Register the reconstructed tracks container
/// @param vTracks Reference to reconstructed tracks container
- void RegisterTracksContainer(L1Vector<CbmL1Track>& vTracks);
+ /// @note If no container is registered, all related routines are omitted
+ void RegisterTracksContainer(L1Vector<CbmL1Track>& vTracks) { fpvTracks = &vTracks; }
/// @brief Sets used detector subsystems
/// @param detID Id of detector
@@ -135,7 +161,11 @@ namespace cbm::ca
template<L1DetectorID DetID>
int ReadHitsForDetector(const TClonesArray* pBrHits);
+ /// @brief Sorts QA hit objects by stations
+ void SortQaHits();
+
/// @brief Saves hit to data structures
+ /// @param hitRecord Filled hit record
///
/// Stores recorded hit information into registered hit containers
void StoreHitRecord(const HitRecord& hitRecord);
@@ -174,19 +204,31 @@ namespace cbm::ca
TClonesArray* fpBrTrdTracks = nullptr; ///< Input branch for reconstructed TRD tracks ("TrdTrack")
TClonesArray* fpBrTofTracks = nullptr; ///< Input branch for reconstructed TOF tracks ("TofTrack")
+
// Pointers to output data containers
L1Vector<CbmL1HitId>* fpvHitIds = nullptr; ///< Pointer to array of hit index objects
- L1Vector<CbmL1HitDebugInfo>* fpvDbgHits = nullptr; ///< Pointer to array of debug hits
+ L1Vector<CbmL1HitDebugInfo>* fpvQaHits = nullptr; ///< Pointer to array of debug hits
L1Vector<CbmL1Track>* fpvTracks = nullptr; ///< Pointer to array of reconstructed tracks
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"};
+ L1Vector<int> vHitMuchIds {"CbmCaTimeSliceReader::vHitMuchIds"};
+ L1Vector<int> vHitTrdIds {"CbmCaTimeSliceReader::vHitTrdIds"};
+ L1Vector<int> vHitTofIds {"CbmCaTimeSliceReader::vHitTofIds"};
// Additional
ECbmTrackingMode fTrackingMode; ///< Tracking mode
- // Indexes
+ // Variables for storing cache
int fNofHits = 0; ///< Stored number of hits
int fNofHitKeys = 0; ///< Recorded number of hit keys
int fFirstHitKey = 0; ///< First index of hit key for the detector subsystem
+
+ std::array<int, L1Constants::size::kMaxNdetectors + 1> fvHitFirstIndexDet = {
+ 0}; ///< First index of hit in detector
};
} // namespace cbm::ca
@@ -206,7 +248,7 @@ int cbm::ca::TimeSliceReader::ReadHitsForDetector(const TClonesArray* pBrHits)
fFirstHitKey = fNofHitKeys;
for (int iH = 0; iH < nHitsTot; ++iH) {
- cbm::ca::HitRecord hitRecord; // Record of hits information
+ HitRecord hitRecord; // Record of hits information
CbmPixelHit* pPixelHit = nullptr; // Pointer to hit object
float phiF = 0.; // Stereo angle of front strips
@@ -269,12 +311,11 @@ int cbm::ca::TimeSliceReader::ReadHitsForDetector(const TClonesArray* pBrHits)
if (ECbmTrackingMode::kMCBM == fTrackingMode && pTofHit->GetZ() > 400) { continue; }
}
- //int iStActive = fpParameters->GetStationIndexActive(iStGeom);
- int iStActive = iStGeom; // !!!! Initialize fParameters
+ int iStActive = fpParameters->GetStationIndexActive(iStGeom, DetID);
if (iStActive == -1) { continue; } // Cut off inactive stations
// Fill out data common for all the detectors
- hitRecord.fStaId = iStGeom;
+ hitRecord.fStaId = iStActive;
hitRecord.fX = pPixelHit->GetX();
hitRecord.fY = pPixelHit->GetY();
hitRecord.fZ = pPixelHit->GetZ();
@@ -287,12 +328,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);
-
- // Update hit external index
- if constexpr (L1DetectorID::kMvd == DetID) { hitRecord.fExtId = -(1 + iH); }
- else {
- hitRecord.fExtId = iH;
- }
+ hitRecord.fExtId = iH;
// Update number of hit keys
if constexpr (L1DetectorID::kSts == DetID) {
diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx
index f54e4cd40c9d28ff74411e2f7b6916eefb191345..c900f001af0ce51c7380ef8321b1955aa271ac8c 100644
--- a/reco/L1/CbmL1.cxx
+++ b/reco/L1/CbmL1.cxx
@@ -1156,8 +1156,11 @@ void CbmL1::DefineSTAPNames(const char* dirName)
fSTAPDataPrefix.ReplaceAll(".root", "");
fSTAPDataPrefix = fSTAPDataPrefix.Strip(TString::EStripType::kBoth, '.');
- TString sDirName = TString(dirName).Length();
- if (!sDirName.Length()) { fSTAPDataDir = pathToRecoOutput.parent_path().string(); }
+ TString sDirName = TString(dirName);
+ if (sDirName.Length() == 0) {
+ fSTAPDataDir = pathToRecoOutput.parent_path().string();
+ L1_SHOW(fSTAPDataDir);
+ }
else if (bfs::exists(sDirName.Data()) && bfs::is_directory(sDirName.Data())) {
fSTAPDataDir = sDirName;
}
diff --git a/reco/L1/CbmL1Hit.h b/reco/L1/CbmL1Hit.h
index 3fabd44e628ced03f5ba365ba294eb174608296b..462548ffce1d66db0ee7ff9b4f2914d87a754692 100644
--- a/reco/L1/CbmL1Hit.h
+++ b/reco/L1/CbmL1Hit.h
@@ -11,9 +11,9 @@
#include "L1Vector.h"
-/// TODO: SZh: Complete the rule of five
-/// TODO: SZh: Make variables private
-/// TODO; SZh: Move class to ca::tools (ca::tools::Hit)
+// TODO: SZh: Complete the rule of five
+// TODO: SZh: Make variables private
+// TODO: SZh: Move class to ca::tools (ca::tools::Hit)
///
/// Identificator for cbm hits with their detector and index in cbm arrays
@@ -23,8 +23,8 @@ public:
CbmL1HitId() = default;
CbmL1HitId(int det, int index) : detId(det), hitId(index) {};
- int detId {0}; ///< detector ID (mvd/sts/etc)
- int hitId {0}; ///< index of hit in the TClonesArray array
+ int detId = -1; ///< detector ID
+ int hitId = -1; ///< index of hit in the TClonesArray array
};
@@ -99,6 +99,7 @@ public:
return msg.str();
}
+ // TODO: SZh 2.03.2023: make the variables private
int ExtIndex; ///< index of hit in the external branch
int iStation; ///< index of station in active stations array
double x; ///< x coordinate of position [cm]
diff --git a/reco/L1/catools/CaToolsMCData.cxx b/reco/L1/catools/CaToolsMCData.cxx
index 15ceb5263cc05b48666fee06918c81103fcc8b46..5a38d5b3354164dee76c7ad2c81fc6bd4fd49de6 100644
--- a/reco/L1/catools/CaToolsMCData.cxx
+++ b/reco/L1/catools/CaToolsMCData.cxx
@@ -26,7 +26,6 @@ MCData::MCData() {}
MCData::MCData(const MCData& other)
: fvPoints(other.fvPoints)
, fvTracks(other.fvTracks)
- , fvPointIndexOfHit(other.fvPointIndexOfHit)
, fmPointLinkMap(other.fmPointLinkMap)
, fmTrackLinkMap(other.fmTrackLinkMap)
{
@@ -61,7 +60,6 @@ void MCData::Swap(MCData& other) noexcept
{
std::swap(fvPoints, other.fvPoints);
std::swap(fvTracks, other.fvTracks);
- std::swap(fvPointIndexOfHit, other.fvPointIndexOfHit);
std::swap(fmPointLinkMap, other.fmPointLinkMap);
std::swap(fmTrackLinkMap, other.fmTrackLinkMap);
}
@@ -88,7 +86,6 @@ void MCData::Clear()
{
fvPoints.clear();
fvTracks.clear();
- fvPointIndexOfHit.clear();
fmPointLinkMap.clear();
fmTrackLinkMap.clear();
}
@@ -126,20 +123,17 @@ std::string MCData::ToString(int verbose) const
if (verbose > 1) {
using std::setfill;
using std::setw;
- const int kNofTracksToPrint = 5;
- msg << "\n Track sample (first 5 tracks):\n";
- msg << setw(10) << setfill(' ') << "ID" << ' ' << setw(10) << setfill(' ') << "motherID" << ' ' << setw(10)
- << setfill(' ') << "pdg" << ' ' << setw(10) << setfill(' ') << "zVTX [cm]" << ' ' << setw(10) << setfill(' ')
- << "px [GeV/c]" << ' ' << setw(10) << setfill(' ') << "py [GeV/c]" << ' ' << setw(10) << setfill(' ')
- << "pz [GeV/c]" << '\n';
+ constexpr int kNofTracksToPrint = 40;
+ constexpr int kNofPointsToPrint = 40;
+ msg << "\n Track sample (first " << kNofTracksToPrint << " tracks):";
+ msg << '\n' << setw(10) << setfill(' ') << fvTracks[0].ToString(verbose, true); // header of the table
for (int i = 0; i < kNofTracksToPrint; ++i) {
- msg << setw(10) << setfill(' ') << fvTracks[i].GetId() << ' ';
- msg << setw(10) << setfill(' ') << fvTracks[i].GetMotherId() << ' ';
- msg << setw(10) << setfill(' ') << fvTracks[i].GetPdgCode() << ' ';
- msg << setw(10) << setfill(' ') << fvTracks[i].GetStartZ() << ' ';
- msg << setw(10) << setfill(' ') << fvTracks[i].GetPx() << ' ';
- msg << setw(10) << setfill(' ') << fvTracks[i].GetPy() << ' ';
- msg << setw(10) << setfill(' ') << fvTracks[i].GetPz() << '\n';
+ msg << '\n' << setw(10) << setfill(' ') << fvTracks[i].ToString(verbose);
+ }
+ msg << "\n Point sample (first " << kNofPointsToPrint << " points):";
+ msg << '\n' << setw(10) << setfill(' ') << fvPoints[0].ToString(verbose, true); // header of the table
+ for (int i = 0; i < kNofPointsToPrint; ++i) {
+ msg << '\n' << setw(10) << setfill(' ') << fvPoints[i].ToString(verbose);
}
}
return msg.str();
diff --git a/reco/L1/catools/CaToolsMCData.h b/reco/L1/catools/CaToolsMCData.h
index 48733204dc2d4c5072f38b18856f513b7097bac3..d0f6ec431d7c8fd01edd104d0865ae8ce81ae2fb 100644
--- a/reco/L1/catools/CaToolsMCData.h
+++ b/reco/L1/catools/CaToolsMCData.h
@@ -107,10 +107,6 @@ namespace ca::tools
/// Gets a reference to the vector of points
const auto& GetPointContainer() const { return fvPoints; }
- /// Gets point index by global index of hit
- /// \param iHit Index of hit
- auto GetPointIndexOfHit(int iHit) const { return fvPointIndexOfHit[iHit]; }
-
/// Gets a reference to MC track by its internal index
const auto& GetTrack(int idx) const { return fvTracks[idx]; }
@@ -130,24 +126,12 @@ namespace ca::tools
/// with hits.
void InitTrackInfo(const L1Vector<CbmL1HitDebugInfo>& vHits);
- /// Registers index of point for a given index of hit
- /// \param iHit Index of hit
- /// \param iPoint Index of point
- void RegisterPointIndexForHit(int iHit, int iPoint)
- {
- assert(int(fvPointIndexOfHit.size()) > iHit);
- fvPointIndexOfHit[iHit] = iPoint;
- }
-
/// Reserves memory for tracks to avoid extra allocations
void ReserveNofTracks(int nTracks) { fvTracks.reserve(nTracks); }
/// Reserves memory for points to avoid extra allocations
void ReserveNofPoints(int nPoints) { fvPoints.reserve(nPoints); }
- /// Reserves total number of used hits in the event
- void ReserveNofHits(int nHits) { fvPointIndexOfHit.reserve(nHits); }
-
/// Prints an example of tracks and points
/// \param verbose Verbose level:
/// - #0: Nothing is printed
@@ -163,9 +147,6 @@ namespace ca::tools
L1Vector<MCPoint> fvPoints = {"ca::tools::MCData::fvPoints"}; ///< Container of points
L1Vector<MCTrack> fvTracks = {"ca::tools::MCData::fvTracks"}; ///< Container of tracks
- /// Correspondence of MC point index to the global hit index
- L1Vector<int> fvPointIndexOfHit = {"ca::tools::MCData::fvPointIndexOfHit"};
-
std::unordered_map<LinkKey, int> fmPointLinkMap = {}; ///< MC point internal index vs. link
std::unordered_map<LinkKey, int> fmTrackLinkMap = {}; ///< MC track internal index vs. link
};
diff --git a/reco/L1/catools/CaToolsMCTrack.cxx b/reco/L1/catools/CaToolsMCTrack.cxx
index 0163c48c92cb9b384a17eabfc5695c1d6c5df33b..4aed3ac5807e61a6112770f8480fffd3a9043263 100644
--- a/reco/L1/catools/CaToolsMCTrack.cxx
+++ b/reco/L1/catools/CaToolsMCTrack.cxx
@@ -11,9 +11,12 @@
#include "CbmL1Hit.h"
+#include <iomanip>
+#include <sstream>
+
#include "CaToolsMCPoint.h"
-using namespace ca::tools;
+using ca::tools::MCTrack;
// ---------------------------------------------------------------------------------------------------------------------
@@ -148,3 +151,69 @@ void MCTrack::SortPointIndexes(const std::function<bool(const int& lhs, const in
{
std::sort(fvPointIndexes.begin(), fvPointIndexes.end(), cmpFn);
}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+std::string MCTrack::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(' ') << fId << ' ';
+ msg << setw(8) << setfill(' ') << fMotherId << ' ';
+ msg << setw(8) << setfill(' ') << fPdgCode << ' ';
+ if (verbose > 1) {
+ msg << setw(8) << setfill(' ') << GetNofHits() << ' ';
+ msg << setw(8) << setfill(' ') << GetNofPoints() << ' ';
+ msg << setw(8) << setfill(' ') << GetNofRecoTracks() << ' ';
+ if (verbose > 2) { msg << setw(8) << setfill(' ') << GetNofTouchTracks() << ' '; }
+ msg << setw(12) << setfill(' ') << GetStartZ() << ' ';
+ msg << setw(12) << setfill(' ') << GetStartT() << ' ';
+ msg << setw(12) << setfill(' ') << GetP() << ' ';
+ }
+ msg << setw(8) << setfill(' ') << IsReconstructable() << ' ';
+ msg << setw(8) << setfill(' ') << IsReconstructed() << ' ';
+ if (verbose > 1) {
+ msg << "\n\t- point indexes: ";
+ for (int index : fvPointIndexes) {
+ msg << index << ' ';
+ }
+ msg << "\n\t- hit indexes: ";
+ for (int index : fvHitIndexes) {
+ msg << index << ' ';
+ }
+ msg << "\n\t- reconstructed track indexes: ";
+ for (int index : fvRecoTrackIndexes) {
+ msg << index << ' ';
+ }
+ if (verbose > 2) {
+ msg << "\n\t- touch track indexes: ";
+ for (int index : fvTouchTrackIndexes) {
+ msg << index << ' ';
+ }
+ }
+ }
+ }
+ }
+ return msg.str();
+}
diff --git a/reco/L1/catools/CaToolsMCTrack.h b/reco/L1/catools/CaToolsMCTrack.h
index 602883d4748061bd01ee657042d37937a9f259a2..338ecaac1a861efa1b3d3844999654f1172c3130 100644
--- a/reco/L1/catools/CaToolsMCTrack.h
+++ b/reco/L1/catools/CaToolsMCTrack.h
@@ -10,6 +10,8 @@
#ifndef CaToolsMCTrack_h
#define CaToolsMCTrack_h 1
+#include "CbmL1Hit.h"
+
#include "TMath.h"
#include <functional>
@@ -310,6 +312,12 @@ namespace ca::tools
/// \param cmpFn Functional object to compare mcPoints
void SortPointIndexes(const std::function<bool(const int& lhs, const int& rhs)>& cmpFn);
+ /// @brief Provides string representation of a track
+ /// @param verbose Verbosity level
+ /// @param header Flag: to print header or data
+ /// @return String representation
+ std::string ToString(int verbose = 1, bool header = false) const;
+
private:
// ****************************
// ** Data variables **
diff --git a/reco/L1/catools/CaToolsQa.cxx b/reco/L1/catools/CaToolsQa.cxx
index 9bd471532625e9e7ca4711a7c473dd4237f7bce0..04fc6bd33012b046d3c25a3842248e18ca543efa 100644
--- a/reco/L1/catools/CaToolsQa.cxx
+++ b/reco/L1/catools/CaToolsQa.cxx
@@ -117,64 +117,65 @@ void Qa::FillHistograms()
// ** Fill distributions for reconstructed tracks (including ghost ones) **
// ************************************************************************
- for (const auto& recoTrk : *fpvRecoTracks) {
- // Reject tracks, which do not contain hits
- if (recoTrk.GetNofHits() < 1) { continue; }
-
- const auto& hitFst = (*fpvHits)[recoTrk.GetHitIndexes()[0]]; // first hit of track
- const auto& hitSnd = (*fpvHits)[recoTrk.GetHitIndexes()[1]]; // second hit of track
- const auto& hitLst = (*fpvHits)[recoTrk.GetHitIndexes()[recoTrk.GetNofHits() - 1]]; // last hit of track
-
- fph_reco_p->Fill(recoTrk.GetP());
- fph_reco_pt->Fill(recoTrk.GetPt());
- fph_reco_phi->Fill(recoTrk.GetPhi());
- fph_reco_nhits->Fill(recoTrk.GetNofHits());
- fph_reco_fsta->Fill(hitFst.GetStationId());
- fph_reco_purity->Fill(100 * recoTrk.GetMaxPurity());
-
- if (recoTrk.GetNDF() > 0) {
- if (recoTrk.IsGhost()) {
- fph_ghost_chi2_ndf->Fill(recoTrk.GetChiSq() / recoTrk.GetNDF());
- fph_ghost_prob->Fill(recoTrk.GetProb());
- }
- else {
- int iTmc = recoTrk.GetMCTrackIndexes()[0]; // Index of first MC-track
- const auto& mcTrk = fpMCData->GetTrack(iTmc);
- // NOTE: reconstructed tracks usually have reference to only one MC-track
- if (mcTrk.IsReconstructable()) {
- fph_reco_chi2_ndf->Fill(recoTrk.GetChiSq() / recoTrk.GetNDF());
- fph_reco_prob->Fill(recoTrk.GetProb());
- }
- else {
- fph_rest_chi2_ndf->Fill(recoTrk.GetChiSq() / recoTrk.GetNDF());
- fph_rest_prob->Fill(recoTrk.GetProb());
- }
- }
- } // 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_nhits->Fill(recoTrk.GetNofHits());
- fph_ghost_fsta->Fill(hitFst.GetStationId());
-
- // 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
- } // loop over recoTrk: end
+ //for (const auto& recoTrk : *fpvRecoTracks) {
+ // // Reject tracks, which do not contain hits
+ // if (recoTrk.GetNofHits() < 1) { continue; }
+
+ // const auto& hitFst = (*fpvHits)[recoTrk.GetHitIndexes()[0]]; // first hit of track
+ // const auto& hitSnd = (*fpvHits)[recoTrk.GetHitIndexes()[1]]; // second hit of track
+ // const auto& hitLst = (*fpvHits)[recoTrk.GetHitIndexes()[recoTrk.GetNofHits() - 1]]; // last hit of track
+
+ // fph_reco_p->Fill(recoTrk.GetP());
+ // fph_reco_pt->Fill(recoTrk.GetPt());
+ // fph_reco_phi->Fill(recoTrk.GetPhi());
+ // fph_reco_nhits->Fill(recoTrk.GetNofHits());
+ // //fph_reco_fsta->Fill(hitFst.GetStationId());
+ // fph_reco_purity->Fill(100 * recoTrk.GetMaxPurity());
+
+ // if (recoTrk.GetNDF() > 0) {
+ // if (recoTrk.IsGhost()) {
+ // fph_ghost_chi2_ndf->Fill(recoTrk.GetChiSq() / recoTrk.GetNDF());
+ // fph_ghost_prob->Fill(recoTrk.GetProb());
+ // }
+ // else {
+ // int iTmc = recoTrk.GetMCTrackIndexes()[0]; // Index of first MC-track
+ // const auto& mcTrk = fpMCData->GetTrack(iTmc);
+ // // NOTE: reconstructed tracks usually have reference to only one MC-track
+ // if (mcTrk.IsReconstructable()) {
+ // fph_reco_chi2_ndf->Fill(recoTrk.GetChiSq() / recoTrk.GetNDF());
+ // fph_reco_prob->Fill(recoTrk.GetProb());
+ // }
+ // else {
+ // fph_rest_chi2_ndf->Fill(recoTrk.GetChiSq() / recoTrk.GetNDF());
+ // fph_rest_prob->Fill(recoTrk.GetProb());
+ // }
+ // }
+ // } // 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_nhits->Fill(recoTrk.GetNofHits());
+ // fph_ghost_fsta->Fill(hitFst.GetStationId());
+ // fph_ghost_fhitR->Fill(hitFst.GetDistFromBP());
+
+ // // 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
+ //} // loop over recoTrk: end
// *************************************
diff --git a/reco/L1/qa/CbmCaOutputQa.cxx b/reco/L1/qa/CbmCaOutputQa.cxx
index b8a7fd83cd164956a031703a617517f60db8563f..51080c550a51ecb2b1c0217e3b91d36d568f38a8 100644
--- a/reco/L1/qa/CbmCaOutputQa.cxx
+++ b/reco/L1/qa/CbmCaOutputQa.cxx
@@ -9,6 +9,7 @@
#include "CbmCaOutputQa.h"
+#include "FairRootManager.h"
#include "Logger.h"
#include "L1InitManager.h"
@@ -27,16 +28,134 @@ void OutputQa::EnableDebugger(const char* filename)
if (!fpDebugger.get()) { fpDebugger = std::make_shared<Debugger>(filename); }
}
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void OutputQa::FillHistograms()
+{
+ // ************************************************************************
+ // ** Fill distributions for reconstructed tracks (including ghost ones) **
+ // ************************************************************************
+
+ for (const auto& recoTrk : fvRecoTracks) {
+ // Reject tracks, which do not contain hits
+ if (recoTrk.GetNofHits() < 1) { continue; }
+
+ const auto& hitFst = (fvHits)[recoTrk.GetHitIndexes()[0]]; // first hit of track
+ const auto& hitSnd = (fvHits)[recoTrk.GetHitIndexes()[1]]; // second hit of track
+ const auto& hitLst = (fvHits)[recoTrk.GetHitIndexes()[recoTrk.GetNofHits() - 1]]; // last hit of track
+
+ // ***************************************************
+ // ** Histograms, requiring reconstructed data only **
+ // ***************************************************
+ //
+ fph_reco_p->Fill(recoTrk.GetP());
+ fph_reco_pt->Fill(recoTrk.GetPt());
+ fph_reco_phi->Fill(recoTrk.GetPhi());
+ fph_reco_tx->Fill(recoTrk.GetTx());
+ fph_reco_ty->Fill(recoTrk.GetTy());
+ fph_reco_nhits->Fill(recoTrk.GetNofHits());
+ fph_reco_fsta->Fill(hitFst.GetStationId());
+ fph_reco_fhitR->Fill(hitFst.GetR());
+
+
+ // ******************************************
+ // ** Histograms, requiring MC information **
+ // ******************************************
+ //
+ if (IsMCUsed()) { fph_reco_purity->Fill(100 * recoTrk.GetMaxPurity()); }
+
+ if (recoTrk.GetNDF() > 0) {
+ if (IsMCUsed() && recoTrk.IsGhost()) {
+ fph_ghost_chi2_ndf->Fill(recoTrk.GetChiSq() / recoTrk.GetNDF());
+ fph_ghost_prob->Fill(recoTrk.GetProb());
+ }
+ else {
+ int iTmc = recoTrk.GetMCTrackIndexes()[0]; // Index of first MC-track
+ const auto& mcTrk = fMCData.GetTrack(iTmc);
+ // NOTE: reconstructed tracks usually have reference to only one MC-track
+ if (mcTrk.IsReconstructable()) {
+ fph_reco_chi2_ndf->Fill(recoTrk.GetChiSq() / recoTrk.GetNDF());
+ fph_reco_prob->Fill(recoTrk.GetProb());
+ }
+ else {
+ fph_rest_chi2_ndf->Fill(recoTrk.GetChiSq() / recoTrk.GetNDF());
+ fph_rest_prob->Fill(recoTrk.GetProb());
+ }
+ }
+ } // 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));
+ }
+
+ 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
+ } // loop over recoTrk: end
+
+
+ // *************************************
+ // ** Fill distributions of MC-tracks **
+ // *************************************
+
+ //for (const auto& mcTrk: fpMCData->GetTrackContainer()) {
+ //
+ //}
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+InitStatus OutputQa::InitCanvases()
+{
+ // ***************************
+ // ** Track distributions
+
+
+ return kSUCCESS;
+}
+
// ---------------------------------------------------------------------------------------------------------------------
//
InitStatus OutputQa::InitDataBranches()
{
+ // Turn off detectors, which hits are not presented in input tree
+ auto* fairManager = FairRootManager::Instance();
+ fbUseMvd = fbUseMvd && fairManager->GetObject("MvdHit");
+ fbUseSts = fbUseSts && fairManager->GetObject("StsHit");
+ fbUseMuch = fbUseMuch && fairManager->GetObject("MuchPixelHit");
+ fbUseTrd = fbUseTrd && fairManager->GetObject("TrdHit");
+ fbUseTof = fbUseTof && fairManager->GetObject("TofHit");
+
+
+ LOG(info) << fName << ": Detector subsystems used:";
+ LOG(info) << "\tMVD: " << (fbUseMvd ? "ON" : "OFF");
+ LOG(info) << "\tSTS: " << (fbUseSts ? "ON" : "OFF");
+ LOG(info) << "\tMuCh: " << (fbUseMuch ? "ON" : "OFF");
+ LOG(info) << "\tTRD: " << (fbUseTrd ? "ON" : "OFF");
+ LOG(info) << "\tTOF: " << (fbUseTof ? "ON" : "OFF");
+
LOG(info) << fName << ": Initializing data branches";
LOG_IF(fatal, !fpParameters.get())
<< fName << ": CA parameters object is not defined. Please, provide initializer or read parameters from binary "
<< "via OutputQa::ReadParameters(filename) from the qa macro";
-
// Initialize IO data manager
if (!fpDataManager.get()) { fpDataManager = std::make_shared<L1IODataManager>(); }
@@ -48,12 +167,12 @@ InitStatus OutputQa::InitDataBranches()
fpTSReader->SetDetector(L1DetectorID::kTrd, fbUseTrd);
fpTSReader->SetDetector(L1DetectorID::kTof, fbUseTof);
- fpTSReader->RegisterIODataManager(fpDataManager);
+ fpTSReader->RegisterParameters(fpParameters);
fpTSReader->RegisterTracksContainer(fvRecoTracks);
- fpTSReader->RegisterHitDebugInfoContainer(fvDbgHits);
- //fpTSReader->RegisterHitIndexContainer(fvHitIds);
+ fpTSReader->RegisterQaHitContainer(fvHits);
+ fpTSReader->RegisterHitIndexContainer(fvHitIds);
- fpTSReader->InitRun();
+ if (!fpTSReader->InitRun()) { return kFATAL; }
// Initialize MC module
if (IsMCUsed()) {
@@ -64,15 +183,70 @@ InitStatus OutputQa::InitDataBranches()
fpMCModule->SetDetector(L1DetectorID::kTrd, fbUseTrd);
fpMCModule->SetDetector(L1DetectorID::kTof, fbUseTof);
- fpMCModule->RegisterInputData(fInputData);
fpMCModule->RegisterMCData(fMCData);
fpMCModule->RegisterRecoTrackContainer(fvRecoTracks);
fpMCModule->RegisterHitIndexContainer(fvHitIds);
+ fpMCModule->RegisterQaHitContainer(fvHits);
fpMCModule->RegisterParameters(fpParameters);
+ fpMCModule->RegisterFirstHitIndexes(fpTSReader->GetHitFirstIndexDet());
- fpMCModule->InitRun();
+ if (!fpMCModule->InitRun()) { return kFATAL; }
}
+ return kSUCCESS;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+InitStatus OutputQa::InitHistograms()
+{
+ fph_reco_p = MakeHistoFromConfig<TH1F>("reco_p");
+ fph_reco_pt = MakeHistoFromConfig<TH1F>("reco_pt");
+ fph_reco_phi = MakeHistoFromConfig<TH1F>("reco_phi");
+ fph_reco_tx = MakeHistoFromConfig<TH1F>("reco_tx");
+ fph_reco_ty = MakeHistoFromConfig<TH1F>("reco_ty");
+ fph_reco_nhits = MakeHistoFromConfig<TH1F>("reco_nhits");
+ fph_reco_fsta = MakeHistoFromConfig<TH1F>("reco_fsta");
+ fph_reco_fhitR = MakeHistoFromConfig<TH1F>("reco_fhitR");
+
+ if (IsMCUsed()) {
+ fph_reco_purity = MakeHistoFromConfig<TH1F>("reco_purity");
+ fph_reco_chi2_ndf = MakeHistoFromConfig<TH1F>("reco_chi2_ndf"); // TODO: Can be filled without MC
+ fph_reco_prob = MakeHistoFromConfig<TH1F>("reco_prob"); // TODO: Can be filled without MC
+ fph_rest_chi2_ndf = MakeHistoFromConfig<TH1F>("rest_chi2_ndf");
+ fph_rest_prob = MakeHistoFromConfig<TH1F>("rest_prob");
+
+ fph_ghost_p = MakeHistoFromConfig<TH1F>("ghost_p");
+ fph_ghost_pt = MakeHistoFromConfig<TH1F>("ghost_pt");
+ fph_ghost_phi = MakeHistoFromConfig<TH1F>("ghost_phi");
+ fph_ghost_nhits = MakeHistoFromConfig<TH1F>("ghost_nhits");
+ fph_ghost_fsta = MakeHistoFromConfig<TH1F>("ghost_fsta");
+ fph_ghost_purity = MakeHistoFromConfig<TH1F>("ghost_purity");
+ fph_ghost_chi2_ndf = MakeHistoFromConfig<TH1F>("ghost_chi2_ndf");
+ fph_ghost_prob = MakeHistoFromConfig<TH1F>("ghost_prob");
+ fph_ghost_tx = MakeHistoFromConfig<TH1F>("ghost_tx");
+ fph_ghost_ty = MakeHistoFromConfig<TH1F>("ghost_ty");
+ fph_ghost_fhitR = MakeHistoFromConfig<TH1F>("ghost_fhitR");
+ fph_ghost_nhits_vs_p = MakeHistoFromConfig<TH2F>("ghost_nhits_vs_p");
+ fph_ghost_fsta_vs_p = MakeHistoFromConfig<TH2F>("ghost_fsta_vs_p");
+ fph_ghost_lsta_vs_fsta = MakeHistoFromConfig<TH2F>("ghost_lsta_vs_fsta");
+ // Residuals and pools of track parameters
+ fph_fst_res_x = MakeHistoFromConfig<TH1F>("fst_res_x");
+ fph_fst_res_y = MakeHistoFromConfig<TH1F>("fst_res_y");
+ fph_fst_res_tx = MakeHistoFromConfig<TH1F>("fst_res_tx");
+ fph_fst_res_ty = MakeHistoFromConfig<TH1F>("fst_res_ty");
+ fph_fst_res_qp = MakeHistoFromConfig<TH1F>("fst_res_qp");
+ fph_fst_res_time = MakeHistoFromConfig<TH1F>("fst_res_time");
+ fph_fst_res_v = MakeHistoFromConfig<TH1F>("fst_res_v");
+
+ fph_fst_pull_x = MakeHistoFromConfig<TH1F>("fst_pull_x");
+ fph_fst_pull_y = MakeHistoFromConfig<TH1F>("fst_pull_y");
+ fph_fst_pull_tx = MakeHistoFromConfig<TH1F>("fst_pull_tx");
+ fph_fst_pull_ty = MakeHistoFromConfig<TH1F>("fst_pull_ty");
+ fph_fst_pull_qp = MakeHistoFromConfig<TH1F>("fst_pull_qp");
+ fph_fst_pull_time = MakeHistoFromConfig<TH1F>("fst_pull_time");
+ fph_fst_pull_v = MakeHistoFromConfig<TH1F>("fst_pull_v");
+ }
return kSUCCESS;
}
@@ -81,30 +255,25 @@ InitStatus OutputQa::InitDataBranches()
//
InitStatus OutputQa::InitTimeSlice()
{
- int nMCTracks = -1;
- int nMCPoints = -1;
- int nHits = -1;
- int nRecoTracks = -1;
-
- // Read hits: fill fvHitIds, fvDbgHits and fpDataManager
- fpTSReader->ReadHits();
- fpDataManager->SendInputData(fInputData);
-
- // Read reconstructed tracks
- fpTSReader->ReadRecoTracks();
+ int nMCTracks = 0;
+ int nMCPoints = 0;
+ int nHits = 0;
+ int nRecoTracks = 0;
- nHits = fInputData.GetNhits();
+ // Read reconstructed input
+ fpTSReader->InitTimeSlice();
+ nHits = fvHits.size();
nRecoTracks = fvRecoTracks.size();
if (IsMCUsed()) {
// Read MC information
fpMCModule->InitEvent(nullptr);
- nMCTracks = fMCData.GetNofTracks();
nMCPoints = fMCData.GetNofPoints();
+ nMCTracks = fMCData.GetNofTracks();
}
- LOG(info) << fName << ": Event or time slice consists from " << nHits << " hits, " << nRecoTracks << " reco tracks, "
- << nMCTracks << " MC tracks, " << nMCPoints << " MC points";
+ LOG_IF(info, fVerbose > 1) << fName << ": Data sample consists of " << nHits << " hits, " << nRecoTracks
+ << " reco tracks, " << nMCTracks << " MC tracks, " << nMCPoints << " MC points";
return kSUCCESS;
}
@@ -119,5 +288,5 @@ void OutputQa::ReadParameters(const char* filename)
manager.ReadParametersObject(filename);
manager.SendParameters(*fpParameters);
- LOG(info) << fpParameters->ToString(10);
+ LOG(info) << fpParameters->ToString(0);
}
diff --git a/reco/L1/qa/CbmCaOutputQa.h b/reco/L1/qa/CbmCaOutputQa.h
index eafee38a1293c0344bc523dd8399ff875b4baa17..91b54aca2da5088aef4321509a01c9060a6a5ceb 100644
--- a/reco/L1/qa/CbmCaOutputQa.h
+++ b/reco/L1/qa/CbmCaOutputQa.h
@@ -13,6 +13,7 @@
#include "CbmCaMCModule.h"
#include "CbmCaTimeSliceReader.h"
#include "CbmL1DetectorID.h"
+#include "CbmL1Hit.h"
#include "CbmQaTask.h"
#include "CaToolsDebugger.h"
@@ -79,20 +80,20 @@ namespace cbm::ca
bool Check() { return true; }
/// @brief Initializes canvases
- InitStatus InitCanvases() { return kSUCCESS; }
+ InitStatus InitCanvases();
/// @brief Initialises data branches in the beginning of the run
InitStatus InitDataBranches();
/// @brief Initializes histograms
- InitStatus InitHistograms() { return kSUCCESS; }
+ InitStatus InitHistograms();
/// @brief Initializes time slice
/// @note Is called in the FairTask::Exec function
InitStatus InitTimeSlice();
- /// @brief Fills histograms
- void FillHistograms() {}
+ /// @brief Fills histograms from time slice or event
+ void FillHistograms();
/// @brief De-initializes histograms
void DeInit() {}
@@ -115,11 +116,60 @@ namespace cbm::ca
std::shared_ptr<::ca::tools::Debugger> fpDebugger = nullptr; ///< Debugger
std::shared_ptr<L1Parameters> fpParameters = nullptr; ///< Tracking parameters object
- L1Vector<CbmL1HitDebugInfo> fvHitIds {"CbmCaOutputQa::fvHitIds"};
- L1Vector<CbmL1HitDebugInfo> fvDbgHits {"CbmCaOutputQa::fvDbgHits"};
+ L1Vector<CbmL1HitId> fvHitIds {"CbmCaOutputQa::fvHitIds"};
+ L1Vector<CbmL1HitDebugInfo> fvHits {"CbmCaOutputQa::fvHits"};
L1Vector<CbmL1Track> fvRecoTracks {"CbmCaOutputQa::fvRecoTracks"};
- L1InputData fInputData; ///< Input hits to CA algo
::ca::tools::MCData fMCData; ///< Input MC data (points and tracks)
+
+
+ // ** Pointers to histogram objects **
+ // 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
+ TH1F* fph_reco_phi = nullptr; ///< Azimuthal 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_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
+ TH1F* fph_rest_prob = nullptr; ///< Fit probability of non-reconstructable tracks
+
+ // Ghost tracks
+ TH1F* fph_ghost_p = nullptr; ///< Total momentum over charge of ghost tracks
+ TH1F* fph_ghost_pt = nullptr; ///< Transverse momentum over charge of ghost tracks
+ TH1F* fph_ghost_phi = nullptr; ///< Azimuthal angle of ghost tracks
+ TH1F* fph_ghost_nhits = nullptr; ///< Hit number of ghost tracks
+ TH1F* fph_ghost_fsta = nullptr; ///< First station index of ghost tracks
+ TH1F* fph_ghost_purity = nullptr; ///< Purity of ghost tracks
+ TH1F* fph_ghost_chi2_ndf = nullptr; ///< Fit chi2 over NDF of ghost tracks
+ TH1F* fph_ghost_prob = nullptr; ///< Fit probability of ghost tracks
+ TH1F* fph_ghost_tx = nullptr; ///< Slope along x-axis of ghost tracks
+ TH1F* fph_ghost_ty = nullptr; ///< Slope along y-axis of ghost tracks
+ TH1F* fph_ghost_fhitR = nullptr; ///< Distance of the first hit from z-axis for ghost tracks
+ TH2F* fph_ghost_nhits_vs_p = nullptr; ///< Hit number vs. total momentum over charge of ghost tracks
+ 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
+
+ // 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]
+ TH1F* fph_fst_res_tx = nullptr; ///< Residual of tx at first track point
+ TH1F* fph_fst_res_ty = nullptr; ///< Residual of ty at first track point
+ TH1F* fph_fst_res_qp = nullptr; ///< Residual of q/p at first track point [GeV/ec]
+ TH1F* fph_fst_res_time = nullptr; ///< Residual of time at first track point [ns]
+ TH1F* fph_fst_res_v = nullptr; ///< Residual of velocity at first track point [c]
+
+ TH1F* fph_fst_pull_x = nullptr; ///< Pull of x at first track point [cm]
+ TH1F* fph_fst_pull_y = nullptr; ///< Pull of y at first track point [cm]
+ TH1F* fph_fst_pull_tx = nullptr; ///< Pull of tx at first track point
+ TH1F* fph_fst_pull_ty = nullptr; ///< Pull of ty at first track point
+ TH1F* fph_fst_pull_qp = nullptr; ///< Pull of q/p at first track point [GeV/ec]
+ TH1F* fph_fst_pull_time = nullptr; ///< Pull of time at first track point [ns]
+ TH1F* fph_fst_pull_v = nullptr; ///< Pull of velocity at first track point [c]
};
} // namespace cbm::ca