From ba851e511dbca8ddf37198ed4668fc9c630190c7 Mon Sep 17 00:00:00 2001
From: "s.zharko@gsi.de" <s.zharko@gsi.de>
Date: Mon, 9 Jan 2023 14:47:36 +0100
Subject: [PATCH] L1: a separate histogrammer for CA tracking
---
reco/L1/CbmCaMCModule.cxx | 39 +++++-
reco/L1/CbmCaMCModule.h | 202 +++++++++++++--------------
reco/L1/CbmL1.cxx | 8 +-
reco/L1/CbmL1.h | 2 +-
reco/L1/CbmL1Hit.h | 94 +++++++++++--
reco/L1/CbmL1Performance.cxx | 13 +-
reco/L1/CbmL1ReadEvent.cxx | 12 +-
reco/L1/CbmL1Track.h | 90 +++++++++---
reco/L1/CbmL1TrackPar.h | 8 +-
reco/L1/catools/CaToolsMCData.cxx | 2 +-
reco/L1/catools/CaToolsMCData.h | 4 +-
reco/L1/catools/CaToolsMCTrack.cxx | 2 +-
reco/L1/catools/CaToolsMCTrack.h | 32 ++---
reco/L1/catools/CaToolsQa.cxx | 212 ++++++++++++++++++++++++++++-
reco/L1/catools/CaToolsQa.h | 186 +++++++++++++++++++++++--
15 files changed, 705 insertions(+), 201 deletions(-)
diff --git a/reco/L1/CbmCaMCModule.cxx b/reco/L1/CbmCaMCModule.cxx
index e76e9e51a8..01b7c99506 100644
--- a/reco/L1/CbmCaMCModule.cxx
+++ b/reco/L1/CbmCaMCModule.cxx
@@ -30,12 +30,13 @@
#include "TLorentzVector.h"
#include "TVector3.h"
+#include <boost/filesystem.hpp>
+
#include <algorithm>
#include <cassert>
#include <fstream> // TODO: SZh 07.12.2022: For debug, should be removed!
#include <stdexcept> // for std::logic_error
-
// *********************************
// ** Action definition functions **
// *********************************
@@ -121,6 +122,23 @@ try {
LOG(info) << "\tTRD: " << (fbUseTrd ? "ON" : "OFF");
LOG(info) << "\tTOF: " << (fbUseTof ? "ON" : "OFF");
+ // Init QA module
+ fpQaModule = std::make_unique<ca::tools::Qa>();
+ fpQaModule->SetParameters(fpParameters);
+ fpQaModule->SetMCData(&fMCData);
+ fpQaModule->SetRecoTrackContainer(fpvRecoTracks);
+ fpQaModule->SetHitContainer(fpvHits);
+
+ // Setup output names
+ // Names are defined using the name for main reconstruction output file, defined from the FairRunAna class
+ boost::filesystem::path sOutPath = FairRunAna::Instance()->GetUserOutputFileName().Data();
+ std::string sOutDir = sOutPath.parent_path().string();
+ if (sOutDir.empty()) { sOutDir = "."; }
+ std::string sOutHistoQa = sOutDir + "/CaHistoQa." + sOutPath.filename().string();
+ fpQaModule->SetOutputFilename(sOutHistoQa.c_str());
+
+ fpQaModule->InitHistograms();
+
return true;
}
catch (const std::logic_error& error) {
@@ -165,7 +183,11 @@ void CbmCaMCModule::InitEvent(CbmEvent* /*pEvent*/)
// ---------------------------------------------------------------------------------------------------------------------
//
-void CbmCaMCModule::InitTrackInfo(const L1Vector<CbmL1Hit>& vHits)
+void CbmCaMCModule::ProcessEvent(CbmEvent*) { fpQaModule->FillHistograms(); }
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void CbmCaMCModule::InitTrackInfo(const L1Vector<CbmL1HitDebugInfo>& vHits)
{
// ----- Initialize stations arrangement and hit indexes
fMCData.InitTrackInfo(vHits);
@@ -208,7 +230,11 @@ void CbmCaMCModule::InitTrackInfo(const L1Vector<CbmL1Hit>& vHits)
// ---------------------------------------------------------------------------------------------------------------------
//
-void CbmCaMCModule::Finish() { std::cout << "\033[1;32mFinishing performance\033[0m\n"; }
+void CbmCaMCModule::Finish()
+{
+ std::cout << "\033[1;32mFinishing performance\033[0m\n";
+ fpQaModule->SaveHistograms();
+}
// **********************************
@@ -217,7 +243,7 @@ void CbmCaMCModule::Finish() { std::cout << "\033[1;32mFinishing performance\033
// ---------------------------------------------------------------------------------------------------------------------
//
-void CbmCaMCModule::MatchPointsWithHits(L1Vector<CbmL1Hit>& vInputExtHits)
+void CbmCaMCModule::MatchPointsWithHits(L1Vector<CbmL1HitDebugInfo>& vInputExtHits)
{
// FIXME: Cleaning up makes it safer, but probably is not needed after old code removal
for (auto& hit : vInputExtHits) {
@@ -383,7 +409,7 @@ int CbmCaMCModule::MatchHitWithMc<L1DetectorID::kTof>(int iHit) const
// ---------------------------------------------------------------------------------------------------------------------
//
-void CbmCaMCModule::MatchRecoAndMCTracks(L1Vector<CbmL1Track>& vRecoTracks, L1Vector<CbmL1Hit>& vInputExtHits)
+void CbmCaMCModule::MatchRecoAndMCTracks(L1Vector<CbmL1Track>& vRecoTracks, L1Vector<CbmL1HitDebugInfo>& vInputExtHits)
{
// NOTE: Potentially there is a case, when
@@ -550,8 +576,7 @@ void CbmCaMCModule::ReadMCPointsForDetector<L1DetectorID::kTof>(CbmMCDataArray*
if (!vbTofPointMatched[iP]) { continue; } // iP was not matched
auto* pExtPoint = L1_DYNAMIC_CAST<CbmTofPoint*>(pPoints->Get(iFile, iEvent, iP));
- LOG_IF(fatal, !pExtPoint) << "NO POINT: TOF: file, event, index = " << iFile << ' ' << iEvent << ' ' << iP
- << '\n';
+ LOG_IF(fatal, !pExtPoint) << "NO POINT: TOF: file, event, index = " << iFile << ' ' << iEvent << ' ' << iP;
if (pExtPoint->GetTrackID() < 0) { continue; } // Point does not have a track
diff --git a/reco/L1/CbmCaMCModule.h b/reco/L1/CbmCaMCModule.h
index 1c7f4ea633..56266d0606 100644
--- a/reco/L1/CbmCaMCModule.h
+++ b/reco/L1/CbmCaMCModule.h
@@ -2,10 +2,10 @@
SPDX-License-Identifier: GPL-3.0-only
Authors: Sergei Zharko [committer] */
-/// \file CbmCaMCModule.h
-/// \brief CA Tracking performance interface for CBM (header)
-/// \since 23.09.2022
-/// \author S.Zharko <s.zharko@gsi.de>
+/// @file CbmCaMCModule.h
+/// @brief CA Tracking performance interface for CBM (header)
+/// @since 23.09.2022
+/// @author S.Zharko <s.zharko@gsi.de>
#ifndef CbmCaMCModule_h
#define CbmCaMCModule_h 1
@@ -41,173 +41,161 @@
class CbmEvent;
class CbmMCDataObject;
-class CbmL1Hit;
+class CbmL1HitDebugInfo;
class CbmL1Track;
enum class L1DetectorID;
-/// Class CbmCaPerformcance is an interface to communicate between
+/// Class CbmCaPerformance is an interface to communicate between
+///
class CbmCaMCModule {
public:
- // *****************************************
- // ** Constructors and destructor **
- // *****************************************
-
- /// Constructor
- /// \param verbosity Verbosity level
+ /// @brief Constructor
+ /// @param verbosity Verbosity level
CbmCaMCModule(int verb = 0, int perfMode = 3) : fVerbose(verb), fPerformanceMode(perfMode) {}
- /// Destructor
+ /// @brief Destructor
~CbmCaMCModule() = default;
- /// Copy constructor
+ /// @brief Copy constructor
CbmCaMCModule(const CbmCaMCModule&) = delete;
- /// Move constructor
+ /// @brief Move constructor
CbmCaMCModule(CbmCaMCModule&&) = delete;
- /// Copy assignment operator
+ /// @brief Copy assignment operator
CbmCaMCModule& operator=(const CbmCaMCModule&) = delete;
- /// Move assignment operator
+ /// @brief Move assignment operator
CbmCaMCModule& operator=(CbmCaMCModule&&) = delete;
-
- // *****************************************
- // ** Action definition functions **
- // *****************************************
-
- /// Defines performance action in the beginning of the run
- /// \return Success flag
- bool InitRun();
-
- /// \brief Defines performance action in the beginning of each event or timeslice
- /// \note This function should be called before hits initialization
- /// \param pEvent Pointer to a current CbmEvent
- void InitEvent(CbmEvent* pEvent);
-
- /// \brief Processes event
- /// Fills histograms and tables, should be called after the tracking done
- void ProcessEvent(CbmEvent* pEvent);
-
- /// \brief Initialize information about arrangement of points and hits of MC tracks within stations and the status
- /// of track reconstructability
- /// \param vHits Vector of hit objects
- /// 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 reconstructablity status
- void InitTrackInfo(const L1Vector<CbmL1Hit>& vHits);
-
- /// Defines performance action in the end of the run
+ /// @brief Defines performance action in the end of the run
void Finish();
- /// Gets the first point index for a given detector subsystem
+ /// @brief Gets the first point index for a given detector subsystem
int GetFirstPointIndex(L1DetectorID detID) const
{
return std::accumulate(fvNofPointsUsed.cbegin(), fvNofPointsUsed.cbegin() + int(detID), 0);
}
- /// Gets the first point index for a given detector subsystem
+ /// @brief Gets the first point index for a given detector subsystem
int GetLastPointIndex(L1DetectorID detID) const
{
return std::accumulate(fvNofPointsUsed.cbegin(), fvNofPointsUsed.cbegin() + int(detID) + 1, 0) - 1;
}
- /// \brief Matches hit with MC point
- /// \tparam DetId Detector ID
- /// \param iHit External index of hit
- /// \return MC-point index in fvMCPoints array
+ /// @brief Gets MC data object
+ /// @return Constant reference to MC data object
+ const ca::tools::MCData& GetMCData() const { return fMCData; }
+
+ /// @brief Gets mutable pointer to MC data object
+ /// @return Mutable reference to MC data object
+ ca::tools::MCData* GetMCData() { return &fMCData; }
+
+ /// @brief Defines performance action in the beginning of each event or time slice
+ /// @note This function should be called before hits initialization
+ /// @param pEvent Pointer to a current CbmEvent
+ void InitEvent(CbmEvent* pEvent);
+
+ /// @brief Defines action on the module in the beginning of the run
+ /// @return Success flag
+ 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);
+
+ /// @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;
- /// Assigns MC point indexes to each hit and hit indexes to each MC point
- /// \param vInputExtHits Vector of external input hits
- void MatchPointsWithHits(L1Vector<CbmL1Hit>& vInputExtHits);
+ /// @brief Assigns MC point indexes to each hit and hit indexes to each MC point
+ /// @param vInputExtHits Vector of external input hits
+ void MatchPointsWithHits(L1Vector<CbmL1HitDebugInfo>& vInputExtHits);
- /// \brief Matches reconstructed tracks with MC tracks
- /// \param vRecoTracks Vector of reconstructed tracks
- /// \param vInputExtHits Vector of external input hits
+ /// @brief Matches reconstructed tracks with MC tracks
+ /// @param vRecoTracks Vector of reconstructed tracks
+ /// @param vInputExtHits Vector of external input hits
+ ///
/// 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(L1Vector<CbmL1Track>& vRecoTracks, L1Vector<CbmL1Hit>& vInputExtHits);
+ void MatchRecoAndMCTracks(L1Vector<CbmL1Track>& vRecoTracks, L1Vector<CbmL1HitDebugInfo>& vInputExtHits);
+ /// @brief Sets used detector subsystems
+ /// @param detID Id of detector
+ /// @param flag Flag: true - detector is used
+ /// @note Should be called before this->Init()
+ void SetDetector(L1DetectorID detID, bool flag);
- // ***********************
- // ** Accessors **
- // ***********************
+ /// @brief Sets pointer to the vector of hits from reference
+ void SetHitContainer(const L1Vector<CbmL1HitDebugInfo>& vHits) { fpvHits = &vHits; }
- /// Gets reference to MC data object
- const ca::tools::MCData& GetMCData() const { return fMCData; }
-
- /// Gets mutable pointer to MC data object
- ca::tools::MCData* GetMCData() { return &fMCData; }
+ /// @brief Sets legacy event mode:
+ /// @param flag Flag:
+ /// - true: runs on events base,
+ /// - false: runs on time-slice base
+ void SetLegacyEventMode(bool flag) { fbLegacyEventMode = flag; }
- /// Registers pointer to the tracking parameters object
+ /// @brief Registers pointer to the tracking parameters object
void SetParameters(const L1Parameters* pParameters) { fpParameters = pParameters; }
- /// Sets used detector subsystems
- /// \note Should be called before this->Init()
- /// \param detID Id of detector
- /// \param flag Flag: true - detector is used
- void SetDetector(L1DetectorID detID, bool flag);
-
- /// Sets legacy event mode:
- /// \param flag: true - runs on events base, false - runs on time-slice base
- void SetLegacyEventMode(bool flag) { fbLegacyEventMode = flag; }
-
- /// Sets pointer to the vector of reconstructed tracks from reference
+ /// @brief Sets pointer to the vector of reconstructed tracks from reference
void SetRecoTrackContainer(const L1Vector<CbmL1Track>& vRecoTracks) { fpvRecoTracks = &vRecoTracks; }
- /// Sets pointer to the vector of hits from reference
- void SetHitContainer(const L1Vector<CbmL1Hit>& vHits) { fpvHits = &vHits; }
-
private:
- // ****************************************
- // ** Local constant expressions **
- // ****************************************
-
- static constexpr int kNdetectors = L1Constants::size::kMaxNdetectors;
-
- // *******************************
- // ** Utility functions **
- // *******************************
-
+ /// @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
- /// \param iPointLocal Local index of MC point
- /// \param detID Detector ID
int CalcGlobMCPointIndex(int iPointLocal, L1DetectorID detID) const
{
return iPointLocal + std::accumulate(fvNofPointsOrig.cbegin(), fvNofPointsOrig.cbegin() + int(detID), 0);
}
- /// Checks class initialization. Throws std::logic_error, if initialization is incomplete at initialization call
+ /// @brief Check class initialization
+ /// @note The function throws std::logic_error, if initialization is incomplete
void CheckInit() const;
- /// Reads MC tracks from external trees and saves them to MCDataObject
+ /// @brief Reads MC tracks from external trees and saves them to MCDataObject
void ReadMCTracks();
- /// Reads MC points from external trees and saves them to MCDataObject
+ /// @brief Reads MC points from external trees and saves them to MCDataObject
void ReadMCPoints();
- /// Reads MC points in particular detector
+ /// @brief Reads MC points in particular detector
template<L1DetectorID DetID>
void ReadMCPointsForDetector(CbmMCDataArray* pPoints);
- /// Fills a single detector-specific MC point
- /// \tparam DetID Detector subsystem ID
- /// \param[in] iExtId Index of point in the external points container
- /// \param[in] iEvent EventID of point in the external points container
- /// \param[in] iFile FileID of point int the external points container
- /// \param[out] intMCPoint Reference to the internal tracking MC point object (ca::tools::MCData)
- /// \return Success flag:
- /// - #true: Point is correct and is to be saved to the MCData object
- /// - #false: Point is incorrect and will be ignored
+ /// @brief Fills a single detector-specific MC point
+ /// @tparam DetID Detector subsystem ID
+ /// @param[in] iExtId Index of point in the external points container
+ /// @param[in] iEvent EventID of point in the external points container
+ /// @param[in] iFile FileID of point int the external points container
+ /// @param[out] intMCPoint Reference to the internal tracking MC point object (ca::tools::MCData)
+ /// @return true Point is correct and is to be saved to the MCData object
+ /// @return false Point is incorrect and will be ignored
template<L1DetectorID DetID>
bool FillMCPoint(int iExtId, int iEvent, int iFile, ca::tools::MCPoint& point);
@@ -248,8 +236,8 @@ private:
CbmMCDataArray* fpTrdPoints = nullptr; ///< TRD MC-points input container
CbmMCDataArray* fpTofPoints = nullptr; ///< TOF MC-points input container
- std::array<int, kNdetectors> fvNofPointsOrig = {}; ///< Number of points by detector provided
- std::array<int, kNdetectors> fvNofPointsUsed = {}; ///< Number of points used in performance
+ 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
// ------ Matches
TClonesArray* fpMvdHitMatches = nullptr; ///< Array of MVD hit matches
@@ -277,7 +265,7 @@ private:
// ********************************************
const L1Vector<CbmL1Track>* fpvRecoTracks = nullptr; ///< Pointer to reconstructed track container
- const L1Vector<CbmL1Hit>* fpvHits = nullptr; ///< Pointer to hits container
+ const L1Vector<CbmL1HitDebugInfo>* fpvHits = nullptr; ///< Pointer to hits container
};
// **********************************************
diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx
index 2a66c06daf..1e533d2c86 100644
--- a/reco/L1/CbmL1.cxx
+++ b/reco/L1/CbmL1.cxx
@@ -978,7 +978,7 @@ void CbmL1::Reconstruct(CbmEvent* event)
t.Tpv[i] = it->Tpv[i];
}
- for (int i = 0; i < L1TrackPar::kNparCov; i++) {
+ for (int i = 0; i < L1TrackPar::kNparCov; ++i) {
t.C[i] = it->CFirst[i];
t.CLast[i] = it->CLast[i];
t.Cpv[i] = it->Cpv[i];
@@ -1009,7 +1009,8 @@ void CbmL1::Reconstruct(CbmEvent* event)
if (fPerformance) {
if (fVerbose > 1) { cout << "Performance..." << endl; }
TrackMatch();
- fpMCModule->MatchRecoAndMCTracks(fvRecoTracks, fvExternalHits);
+ //fpMCModule->MatchRecoAndMCTracks(fvRecoTracks, fvHitDebugInfo);
+ //fpMCModule->ProcessEvent(event);
EfficienciesPerformance();
HistoPerformance();
TrackFitPerformance();
@@ -1024,9 +1025,12 @@ void CbmL1::Reconstruct(CbmEvent* event)
// ----- Finish CbmStsFitPerformanceTask task -----------------------------
void CbmL1::Finish()
{
+ //if (fPerformance) { fpMCModule->Finish(); }
+
TDirectory* curr = gDirectory;
TFile* currentFile = gFile;
+
// Open output file and write histograms
boost::filesystem::path p = (FairRunAna::Instance()->GetUserOutputFileName()).Data();
std::string dir = p.parent_path().string();
diff --git a/reco/L1/CbmL1.h b/reco/L1/CbmL1.h
index 5cd6d31c20..2a45101085 100644
--- a/reco/L1/CbmL1.h
+++ b/reco/L1/CbmL1.h
@@ -449,7 +449,7 @@ private:
L1InitManager fInitManager; ///< Tracking parameters data manager
L1IODataManager fIODataManager; ///< Input-output data manager
- std::unique_ptr<CbmCaMCModule> fpMCModule = nullptr; ///< MC-module for tracking
+ //std::unique_ptr<CbmCaMCModule> fpMCModule = nullptr; ///< MC-module for tracking
public:
diff --git a/reco/L1/CbmL1Hit.h b/reco/L1/CbmL1Hit.h
index 02007d3c92..3fabd44e62 100644
--- a/reco/L1/CbmL1Hit.h
+++ b/reco/L1/CbmL1Hit.h
@@ -5,8 +5,16 @@
#ifndef _CbmL1Hit_h_
#define _CbmL1Hit_h_
+#include <iomanip>
+#include <sstream>
+#include <string>
+
#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)
+
///
/// Identificator for cbm hits with their detector and index in cbm arrays
///
@@ -25,16 +33,82 @@ public:
///
class CbmL1HitDebugInfo { // TODO: SZh 21.09.2022: Replace instances of this class with L1Hit
public:
- 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]
- double y; ///< y coordinate of position [cm]
- double time; ///< hit time [ns]
- double dx; ///< x coordinate error [cm]
- double dy; ///< y coordinate error [cm]
- double dt; ///< time error [ns]
- double dxy; ///< covariance between x and y [cm2]
- int Det; ///< detector subsystem ID
+ /// Gets detector type
+ /// 0 - MVD
+ /// 1 - STS
+ /// 2 - MuCh
+ /// 3 - TRD
+ /// 4 - TOF
+ // TODO: SZh 02.03.2023: Replace det. ID with L1DetectorID
+ int GetDetectorType() const { return Det; }
+
+ /// Gets distance from z-axis [cm]
+ float GetR() const { return std::sqrt(x * x + y * y); }
+
+ /// Gets index of the hit in the external container
+ int GetExternalHitId() const { return ExtIndex; }
+
+ /// Gets MC point index container
+ const auto& GetMCPointIndexes() const { return mcPointIds; }
+
+ /// Gets index of matched MC point
+ int GetMCPointIndex() const { return mcPointIds.size() ? mcPointIds[0] : -1; }
+
+ /// Gets global index of active tracking station
+ int GetStationId() const { return iStation; }
+
+ /// Gets time measurement [ns]
+ float GetT() const { return time; }
+
+ /// Gets x component of position [cm]
+ float GetX() const { return x; }
+
+ /// Gets y component of position [cm]
+ float GetY() const { return y; }
+
+ /// @brief String representation of the object
+ /// @param verbose Verbosity level
+ /// @param header If true, header will be printed
+ std::string ToString(int verbose = 10, bool header = false) const
+ {
+ using std::setfill;
+ using std::setw;
+ std::stringstream msg;
+ if (header) {
+ msg << setw(8) << setfill(' ') << "ext. ID" << ' ';
+ msg << setw(8) << setfill(' ') << "st. ID" << ' ';
+ msg << setw(12) << setfill(' ') << "x [cm]" << ' ';
+ msg << setw(12) << setfill(' ') << "y [cm]" << ' ';
+ msg << setw(12) << setfill(' ') << "time [ns]" << ' ';
+ msg << setw(8) << setfill(' ') << "Det. ID" << ' ';
+ }
+ else {
+ msg << setw(8) << setfill(' ') << ExtIndex << ' ';
+ msg << setw(8) << setfill(' ') << iStation << ' ';
+ msg << setw(12) << setfill(' ') << x << ' ';
+ msg << setw(12) << setfill(' ') << y << ' ';
+ msg << setw(12) << setfill(' ') << time << ' ';
+ msg << setw(8) << setfill(' ') << Det << ' ';
+ if (verbose > 0) {
+ msg << "\n\tMC point indexes: ";
+ for (int iP : mcPointIds) {
+ msg << iP << ' ';
+ }
+ }
+ }
+ return msg.str();
+ }
+
+ 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]
+ double y; ///< y coordinate of position [cm]
+ double time; ///< hit time [ns]
+ double dx; ///< x coordinate error [cm]
+ double dy; ///< y coordinate error [cm]
+ double dt; ///< time error [ns]
+ double dxy; ///< covariance between x and y [cm2]
+ int Det; ///< detector subsystem ID
L1Vector<int> mcPointIds {"CbmL1HitDebugInfo::mcPointIds"}; ///< indices of CbmL1MCPoint in L1->fvMCPoints array
};
diff --git a/reco/L1/CbmL1Performance.cxx b/reco/L1/CbmL1Performance.cxx
index 38432f6307..a7b56af548 100644
--- a/reco/L1/CbmL1Performance.cxx
+++ b/reco/L1/CbmL1Performance.cxx
@@ -265,14 +265,16 @@ struct TL1PerfEfficiencies : public TL1Efficiencies {
{
L1_assert(nEvents != 0);
int NCounters = mc.GetNcounters();
- std::vector<std::string> rowNames(20);
+ std::vector<std::string> rowNames(NCounters + 2);
for (int iC = 0; iC < NCounters; ++iC) {
rowNames[iC] = std::string(names[iC].Data());
}
+ rowNames[NCounters] = "Ghost prob.";
+ rowNames[NCounters + 1] = "N ghosts";
std::vector<std::string> colNames = {"Eff.", "Killed", "Length", "Fakes", "Clones",
"All Reco", "All MC", "MCl(hits)", "MCl(MCps)"};
- CbmQaTable* aTable = new CbmQaTable(nameOfTable.c_str(), "Track Efficiency", 20, 9);
+ CbmQaTable* aTable = new CbmQaTable(nameOfTable.c_str(), "Track Efficiency", NCounters + 2, 9);
aTable->SetNamesOfRows(rowNames);
aTable->SetNamesOfCols(colNames);
for (int iC = 0; iC < NCounters; iC++) {
@@ -286,6 +288,9 @@ struct TL1PerfEfficiencies : public TL1Efficiencies {
aTable->SetCell(iC, 7, mc_length_hits.counters[iC] / double(mc.counters[iC]));
aTable->SetCell(iC, 8, mc_length.counters[iC] / double(mc.counters[iC]));
}
+ aTable->SetCell(NCounters, 0, ratio_ghosts);
+ aTable->SetCell(NCounters + 1, 0, ghosts);
+
if (ifPrintTableToLog) {
cout << *aTable; // print a table to log
}
@@ -293,8 +298,6 @@ struct TL1PerfEfficiencies : public TL1Efficiencies {
else {
delete aTable;
}
-
- cout << "Ghost probability : " << ratio_ghosts << " | " << ghosts << endl;
};
TL1TracksCatCounters<double> ratio_killed;
@@ -386,7 +389,7 @@ void CbmL1::EfficienciesPerformance()
int mc_length = mtra.NMCStations();
if (reco) {
for (unsigned int irt = 0; irt < rTracks.size(); irt++) {
- ratio_length += static_cast<double>(rTracks[irt]->GetNOfHits()) * rTracks[irt]->GetMaxPurity() / mc_length_hits;
+ ratio_length += static_cast<double>(rTracks[irt]->GetNofHits()) * rTracks[irt]->GetMaxPurity() / mc_length_hits;
ratio_fakes += 1 - rTracks[irt]->GetMaxPurity();
}
}
diff --git a/reco/L1/CbmL1ReadEvent.cxx b/reco/L1/CbmL1ReadEvent.cxx
index 5b7b083f9c..4925155135 100644
--- a/reco/L1/CbmL1ReadEvent.cxx
+++ b/reco/L1/CbmL1ReadEvent.cxx
@@ -1153,7 +1153,7 @@ void CbmL1::ReadEvent(CbmEvent* event)
fIODataManager.ResetInputData();
fIODataManager.ReserveNhits(nHits);
fIODataManager.SetNhitKeys(NStrips);
- if (fPerformance) { fpMCModule->GetMCData()->ReserveNofHits(nHits); }
+ //if (fPerformance) { fpMCModule->GetMCData()->ReserveNofHits(nHits); }
// ----- Fill
for (int iHit = 0; iHit < nHits; ++iHit) {
@@ -1196,10 +1196,10 @@ void CbmL1::ReadEvent(CbmEvent* event)
fvHitDebugInfo.push_back(s);
fvHitPointIndexes.push_back(th.iMC);
- fpMCModule->GetMCData()->RegisterPointIndexForHit(iHit, th.iMC);
+ //if (fPerformance) { fpMCModule->GetMCData()->RegisterPointIndexForHit(iHit, th.iMC); }
}
- if (fPerformance) { HitMatch(); } /// OLD
- if (fPerformance) { fpMCModule->MatchPointsWithHits(fvExternalHits); } /// NEW
+ if (fPerformance) { HitMatch(); } /// OLD
+ //if (fPerformance) { fpMCModule->MatchPointsWithHits(fvHitDebugInfo); } /// NEW
if (fVerbose >= 2) cout << "ReadEvent: mvd and sts are saved." << endl;
@@ -1468,7 +1468,7 @@ bool CbmL1::ReadMCPoint(CbmL1MCPoint* MC, int iPoint, int file, int event, int i
void CbmL1::HitMatch()
{
// Clear contents
- for (auto& hit : fvExternalHits) {
+ for (auto& hit : fvHitDebugInfo) {
hit.mcPointIds.clear();
}
@@ -1477,7 +1477,7 @@ void CbmL1::HitMatch()
}
// Fill new contents
- const int NHits = fvExternalHits.size();
+ const int NHits = fvHitDebugInfo.size();
for (int iH = 0; iH < NHits; iH++) {
CbmL1HitDebugInfo& hit = fvHitDebugInfo[iH];
int iP = fvHitPointIndexes[iH];
diff --git a/reco/L1/CbmL1Track.h b/reco/L1/CbmL1Track.h
index b7434e01bc..6ea035c7ed 100644
--- a/reco/L1/CbmL1Track.h
+++ b/reco/L1/CbmL1Track.h
@@ -1,6 +1,6 @@
-/* Copyright (C) 2006-2021 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+/* Copyright (C) 2006-2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
SPDX-License-Identifier: GPL-3.0-only
- Authors: Ivan Kisel, Sergey Gorbunov, Denis Bertini [committer], Igor Kulakov, Maksym Zyzak */
+ Authors: Ivan Kisel, Sergey Gorbunov, Denis Bertini [committer], Igor Kulakov, Maksym Zyzak, Sergei Zharko */
/*
*====================================================================
@@ -26,6 +26,8 @@
#include "CbmL1MCTrack.h"
#include "CbmL1TrackPar.h"
+#include "TMath.h"
+
#include <map>
#include <vector>
using std::map;
@@ -37,8 +39,7 @@ class CbmL1Track : public CbmL1TrackPar {
public:
CbmL1Track() : Hits(), nStations(0), index(0), fTrackTime(0.), hitMap(), mcTracks(), maxPurity(-1) {}
- int GetNOfHits() { return Hits.size(); }
-
+ /// Adds pointer to MC track (TODO: remove this and related methods)
void AddMCTrack(CbmL1MCTrack* mcTr) { mcTracks.push_back(mcTr); }
/// Adds an index of MC track index
@@ -51,17 +52,69 @@ public:
fvMcTrackIndexes.clear();
}
- vector<CbmL1MCTrack*>& GetMCTracks() { return mcTracks; }
- CbmL1MCTrack* GetMCTrack() { return mcTracks[0]; }
- int GetNMCTracks() { return mcTracks.size(); }
- bool IsGhost() { return (maxPurity < CbmL1Constants::MinPurity); }
+ /// Gets Chi-square of track fit model
+ double GetChiSq() const { return chi2; }
+
+ /// Gets a reference to hit indexes array
+ const vector<int>& GetHitIndexes() const { return Hits; }
+
+ bool IsGhost() const { return (maxPurity < CbmL1Constants::MinPurity); }
+
+ /// Gets NDF of track fit model
+ int GetNDF() const { return NDF; }
+ /// Gets number of MC tracks
+ // TODO: Remove this method
+ int GetNMCTracks() const { return mcTracks.size(); }
/// Gets a reference to MC track indexes
const auto& GetMCTrackIndexes() const { return fvMcTrackIndexes; }
/// Gets max purity
- double GetMaxPurity() { return maxPurity; }
+ double GetMaxPurity() const { return maxPurity; }
+
+ /// Gets container of pointers to MC tracks
+ // TODO: this function is to be replaced with GetMCTrackIndexes()
+ vector<CbmL1MCTrack*>& GetMCTracks() { return mcTracks; }
+
+ /// Gets pointer of the associated MC track
+ CbmL1MCTrack* GetMCTrack() { return mcTracks[0]; }
+
+ /// Gets number of hits of the track
+ int GetNofHits() const { return Hits.size(); }
+
+ /// Gets number of associated MC tracks
+ int GetNofMCTracks() const { return fvMcTrackIndexes.size(); }
+
+ /// Gets number of stations
+ int GetNofStations() const { return nStations; }
+
+ /// Gets absolute value of momentum divided by the charge (absolute value) of particle [GeV/ec]
+ double GetP() const { return fabs(T[4]) > 1.e-10 ? 1. / fabs(T[4]) : 1.e-10; }
+
+ /// Gets transverse momentum
+ double GetPt() const { return sqrt(GetP() * GetP() - GetPz() * GetPz()); }
+
+ /// Gets azimuthal angle of the track
+ double GetPhi() const { return TMath::ATan2(-GetTy(), -GetTx()); }
+
+ /// Gets probability of track fit model
+ double GetProb() const { return TMath::Prob(chi2, NDF); }
+
+ /// Gets x-component of momentum divided by the charge (absolute value) of particle [GeV/ec]
+ double GetPx() const { return GetPz() * GetTx(); }
+
+ /// Gets y-component of momentum divided by the charge (absolute value) of particle [GeV/ec]
+ double GetPy() const { return GetPz() * GetTy(); }
+
+ /// Gets z-component of momentum divided by the charge (absolute value) of particle [GeV/ec]
+ double GetPz() const { return std::sqrt(GetP() * GetP() / (GetTx() * GetTx() + GetTy() * GetTy() + 1)); }
+
+ /// Gets track slope along x-axis
+ double GetTx() const { return T[2]; }
+
+ /// Gets track slope along y-axis
+ double GetTy() const { return T[3]; }
/// Sets max purity
/// NOTE: max purity is calculated as a ratio of max number of hits left by an actual track and the
@@ -70,21 +123,26 @@ public:
static bool compareChi2(const CbmL1Track& a, const CbmL1Track& b) { return (a.chi2 < b.chi2); }
+
static bool comparePChi2(const CbmL1Track* a, const CbmL1Track* b) { return (a->chi2 < b->chi2); }
- double Tpv[L1TrackPar::kNparTr], Cpv[L1TrackPar::kNparCov];
+ double Tpv[L1TrackPar::kNparTr]; ///< Track parameters at primary vertex
+ double Cpv[L1TrackPar::kNparCov]; ///< Track covariance matrix at primary vertex
+ double TLast[L1TrackPar::kNparTr]; ///< Track parameters in the end of the track
+ double CLast[L1TrackPar::kNparCov]; ///< Track covariance matrix at the end of the track
- double TLast[L1TrackPar::kNparTr], CLast[L1TrackPar::kNparCov];
- vector<int> Hits;
- int nStations;
- int index;
+ vector<int> Hits; ///< Indexes of hits of this track
+ int nStations; ///< Number of stations with hits of this track
+ int index; ///< Index of this track (TODO: it seems to be not initialized)
- double fTrackTime;
+ double fTrackTime; ///< Time of the track [ns] ???
map<int, int> hitMap; // N hits (second) from each mcTrack (first is a MC track ID) belong to current recoTrack
// FIXME: SZh 14.12.2022: map => unordered_map
+
+
private:
// next members filled and used in Performance
vector<CbmL1MCTrack*> mcTracks; // array of assosiated recoTracks. Should be only one.
@@ -92,7 +150,7 @@ private:
L1Vector<int> fvMcTrackIndexes = {"CbmL1Track::fvMcTrackIndexes"}; // global indexes of MC tracks
// NOTE: mcTracks should be replaced with fvMcTrackIndexes
- double maxPurity; // maximum persent of hits, which belong to one mcTrack.
+ double maxPurity; ///< Maximum persent of hits, which belong to one mcTrack.
};
#endif
diff --git a/reco/L1/CbmL1TrackPar.h b/reco/L1/CbmL1TrackPar.h
index 57fe79607b..421d53e03c 100644
--- a/reco/L1/CbmL1TrackPar.h
+++ b/reco/L1/CbmL1TrackPar.h
@@ -16,10 +16,12 @@ public:
double* GetCovMatrix() { return C; }
double& GetRefChi2() { return chi2; }
int& GetRefNDF() { return NDF; }
- double GetMass() { return mass; }
- bool IsElectron() { return is_electron; }
+ double GetMass() const { return mass; }
+ bool IsElectron() const { return is_electron; }
- double T[L1TrackPar::kNparTr], C[L1TrackPar::kNparCov], chi2;
+ double T[L1TrackPar::kNparTr];
+ double C[L1TrackPar::kNparCov];
+ double chi2;
int NDF;
double mass; // mass hypothesis
bool is_electron;
diff --git a/reco/L1/catools/CaToolsMCData.cxx b/reco/L1/catools/CaToolsMCData.cxx
index 64573973b6..c01fe96611 100644
--- a/reco/L1/catools/CaToolsMCData.cxx
+++ b/reco/L1/catools/CaToolsMCData.cxx
@@ -95,7 +95,7 @@ void MCData::Clear()
// ---------------------------------------------------------------------------------------------------------------------
//
-void MCData::InitTrackInfo(const L1Vector<CbmL1Hit>& vHits)
+void MCData::InitTrackInfo(const L1Vector<CbmL1HitDebugInfo>& vHits)
{
for (auto& aTrk : fvTracks) {
// Assign hits to tracks
diff --git a/reco/L1/catools/CaToolsMCData.h b/reco/L1/catools/CaToolsMCData.h
index 1fab4a3125..48733204dc 100644
--- a/reco/L1/catools/CaToolsMCData.h
+++ b/reco/L1/catools/CaToolsMCData.h
@@ -128,7 +128,7 @@ namespace ca::tools
/// Initialize tracks: defines indexes of hits and points related to the track, 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(const L1Vector<CbmL1Hit>& vHits);
+ void InitTrackInfo(const L1Vector<CbmL1HitDebugInfo>& vHits);
/// Registers index of point for a given index of hit
/// \param iHit Index of hit
@@ -146,7 +146,7 @@ namespace ca::tools
void ReserveNofPoints(int nPoints) { fvPoints.reserve(nPoints); }
/// Reserves total number of used hits in the event
- void ReserveNofHits(int nHits) { fvPointIndexOfHit.reset(nHits); }
+ void ReserveNofHits(int nHits) { fvPointIndexOfHit.reserve(nHits); }
/// Prints an example of tracks and points
/// \param verbose Verbose level:
diff --git a/reco/L1/catools/CaToolsMCTrack.cxx b/reco/L1/catools/CaToolsMCTrack.cxx
index 1a07b81dee..0163c48c92 100644
--- a/reco/L1/catools/CaToolsMCTrack.cxx
+++ b/reco/L1/catools/CaToolsMCTrack.cxx
@@ -29,7 +29,7 @@ void MCTrack::Clear()
// ---------------------------------------------------------------------------------------------------------------------
//
-void MCTrack::InitHitsInfo(const L1Vector<CbmL1Hit>& vHits)
+void MCTrack::InitHitsInfo(const L1Vector<CbmL1HitDebugInfo>& vHits)
{
// NOTE: vHits must be sorted over stations!
fMaxNofHitsOnStation = 0;
diff --git a/reco/L1/catools/CaToolsMCTrack.h b/reco/L1/catools/CaToolsMCTrack.h
index da2a1bcef8..602883d474 100644
--- a/reco/L1/catools/CaToolsMCTrack.h
+++ b/reco/L1/catools/CaToolsMCTrack.h
@@ -10,13 +10,16 @@
#ifndef CaToolsMCTrack_h
#define CaToolsMCTrack_h 1
+#include "TMath.h"
+
#include <functional>
#include "CaToolsLinkKey.h"
+#include "L1Constants.h"
#include "L1Undef.h"
#include "L1Vector.h"
-class CbmL1Hit;
+class CbmL1HitDebugInfo;
namespace ca::tools
{
@@ -142,7 +145,7 @@ namespace ca::tools
int GetPdgCode() const { return fPdgCode; }
/// Gets azimuthal angle [rad]
- double GetPhi() const;
+ double GetPhi() const { return TMath::ATan2(-fMom[1], -fMom[0]); }
/// Gets a reference to associated point indexes
const auto& GetPointIndexes() const { return fvPointIndexes; }
@@ -189,6 +192,12 @@ namespace ca::tools
/// Gets total number of stations with MC points
int GetTotNofStationsWithPoint() const { return fTotNofStationsWithPoint; }
+ /// Gets track slope along x-axis
+ double GetTx() const { return fMom[0] / fMom[2]; }
+
+ /// Gets track slope along y-axis
+ double GetTy() const { return fMom[1] / fMom[2]; }
+
/// Gets a reference to vector of reconstructed track indexes, not associated with this MC track but containing some
/// hits, produced by this MC track
const auto& GetTouchTrackIndexes() const { return fvTouchTrackIndexes; }
@@ -199,7 +208,7 @@ namespace ca::tools
/// #2) Maximal number of hits within one station
/// #3) Number of consecutive stations with a hit in MC track
/// \param vHits Vector of hits for a given TS
- void InitHitsInfo(const L1Vector<CbmL1Hit>& vHits);
+ void InitHitsInfo(const L1Vector<CbmL1HitDebugInfo>& vHits);
/// \brief Initializes information about MC track points arrangement within stations
/// Defines:
@@ -345,21 +354,4 @@ namespace ca::tools
} // namespace ca::tools
-
-// ************************************
-// ** Inline implementations **
-// ************************************
-
-inline double ca::tools::MCTrack::GetPhi() const
-{
- // NOTE: Implementation was taken from ROOT TMath::Atan2
- constexpr double kRootPi = 3.14159265358979323846; // Value of Pi, used in ROOT TMath
- if (fMom[0] != 0) { return std::atan2(-fMom[1], -fMom[0]); }
- if (fMom[1] == 0) { return 0; }
- if (fMom[1] < 0) { return kRootPi / 2; }
- else {
- return -kRootPi / 2;
- }
-}
-
#endif // CaToolsMCTrack_h
diff --git a/reco/L1/catools/CaToolsQa.cxx b/reco/L1/catools/CaToolsQa.cxx
index 1f693a8995..9bd4715326 100644
--- a/reco/L1/catools/CaToolsQa.cxx
+++ b/reco/L1/catools/CaToolsQa.cxx
@@ -1,6 +1,6 @@
-/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+/* Copyright (C) 2022-2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
SPDX-License-Identifier: GPL-3.0-only
- Authors: Sergey Gorbunov, Sergei Zharko [committer] */
+ Authors: Sergei Zharko [committer] */
/// \file CaToolsQa.h
/// \brief Tracking performance class (implementation)
@@ -9,18 +9,220 @@
#include "CaToolsQa.h"
+#include "CbmL1Hit.h" // TODO: SZh: move CbmL1Hit to ca::tools::Hit
+#include "CbmL1Track.h" // TODO: SZh: move CbmL1Track to ca::tools::RecoTrack
+
+#include "FairRunAna.h"
+
#include <utility> // for std::move
+#include "CaToolsMCData.h"
+#include "L1Parameters.h"
+
using namespace ca::tools;
// ---------------------------------------------------------------------------------------------------------------------
//
-void Qa::InitHistograms() {}
+void Qa::InitHistograms()
+{
+ LOG_IF(fatal, !fpOutFile) << "CA QA: Output file was not defined";
+
+ TDirectory* currentDir = gDirectory;
+ gDirectory = fpHistoDir;
+
+ // ----- Initialize histograms and profiles --------------------------------------------------------------------------
+ // Distributions of reconstructed tracks by different quantities
+ fph_reco_p = MakeHisto<TH1F>("ca_reco_p", "", kNbinsRecoP, 0., kAxMaxRecoP);
+ fph_reco_pt = MakeHisto<TH1F>("ca_reco_pt", "", kNbinsRecoPt, 0., kAxMaxRecoPt);
+ fph_reco_phi = MakeHisto<TH1F>("ca_reco_phi", "", kNbinsRecoPhi, -kAxMaxRecoPhi, kAxMaxRecoPhi);
+ fph_reco_nhits = MakeHisto<TH1F>("ca_reco_nhits", "", kNbinsNofHits, -0.5, kNbinsNofHits - 0.5);
+ fph_reco_fsta = MakeHisto<TH1F>("ca_reco_fsta", "", kNbinsNofSta, -0.5, kNbinsNofSta - 0.5);
+ fph_reco_purity = MakeHisto<TH1F>("ca_reco_purity", "", kNbinsPurity, -0.5, 100.5);
+ fph_reco_chi2_ndf = MakeHisto<TH1F>("ca_reco_chi2_ndf", "", kNbinsChi2NDF, -0.5, kAxMaxChi2NDF);
+ fph_reco_prob = MakeHisto<TH1F>("ca_reco_prob", "", kNbinsProb, 0., 1.01);
+ fph_rest_chi2_ndf = MakeHisto<TH1F>("ca_rest_chi2_ndf", "", kNbinsChi2NDF, -0.5, kAxMaxChi2NDF);
+ fph_rest_prob = MakeHisto<TH1F>("ca_rest_prob", "", kNbinsProb, 0., 1.01);
+
+ fph_reco_p->SetTitle("Momentum of reconstructed track;(p/q)_{reco} [GeV/ec]");
+ fph_reco_pt->SetTitle("Transverse momentum of reconstructed track;(p_{T}/q)_{reco} [GeV/ec]");
+ fph_reco_phi->SetTitle("Azimuthal angle of reconstructed track;#phi_{reco} [rad]");
+ fph_reco_nhits->SetTitle("Number of hits in a reconstructed track;N_{hits}");
+ fph_reco_fsta->SetTitle("First station of reconstructed track;ID_{f.st.}");
+ fph_reco_purity->SetTitle("Purity of reconstructed track;purity [%]");
+ fph_reco_chi2_ndf->SetTitle("Fit #chi^{2}/NDF of reconstructed track;(#chi^{2}/NDF)_{reco}");
+ fph_reco_chi2_ndf->SetTitle("Fit #chi^{2}/NDF of the rest track;(#chi^{2}/NDF)_{reco}");
+ fph_reco_prob->SetTitle("Fit probability of reconstructed track;Prob_{reco}");
+ fph_rest_prob->SetTitle("Fit probability of the rest tracks;Prob_{reco}");
+
+
+ // Distributions of reconstructed tracks by different quantities
+ fph_ghost_p = MakeHisto<TH1F>("ca_ghost_p", "", kNbinsRecoP, 0., kAxMaxRecoP);
+ fph_ghost_pt = MakeHisto<TH1F>("ca_ghost_pt", "", kNbinsRecoPt, 0., kAxMaxRecoPt);
+ fph_ghost_phi = MakeHisto<TH1F>("ca_ghost_phi", "", kNbinsRecoPhi, -kAxMaxRecoPhi, kAxMaxRecoPhi);
+ fph_ghost_nhits = MakeHisto<TH1F>("ca_ghost_nhits", "", kNbinsNofHits, -0.5, kNbinsNofHits - 0.5);
+ fph_ghost_fsta = MakeHisto<TH1F>("ca_ghost_fsta", "", kNbinsNofSta, -0.5, kNbinsNofSta - 0.5);
+ fph_ghost_purity = MakeHisto<TH1F>("ca_ghost_purity", "", kNbinsPurity, -0.5, 100.5);
+ fph_ghost_chi2_ndf = MakeHisto<TH1F>("ca_ghost_chi2_ndf", "", kNbinsChi2NDF, 0., kAxMaxChi2NDF);
+ fph_ghost_prob = MakeHisto<TH1F>("ca_ghost_prob", "", kNbinsProb, 0., 1.01);
+ fph_ghost_tx = MakeHisto<TH1F>("ca_ghost_tx", "", kNbinsTx, -2., 2.);
+ fph_ghost_ty = MakeHisto<TH1F>("ca_ghost_ty", "", kNbinsTy, -2., 2.);
+ fph_ghost_fhitR = MakeHisto<TH1F>("ca_ghost_fhitR", "", kNbinsHitDist, 0., 150.);
+ fph_ghost_nhits_vs_p =
+ MakeHisto<TH2F>("ca_ghost_nhits_vs_p", "", kNbinsRecoP, 0., kAxMaxRecoP, kNbinsNofHits, -0.5, kNbinsNofHits - 0.5);
+ fph_ghost_fsta_vs_p =
+ MakeHisto<TH2F>("ca_ghost_fsta_vs_p", "", kNbinsRecoP, 0., kAxMaxRecoP, kNbinsNofSta, -0.5, kNbinsNofSta - 0.5);
+ fph_ghost_lsta_vs_fsta = MakeHisto<TH2F>("ca_ghost_lsta_vs_fsta", "", kNbinsNofSta, -0.5, kNbinsNofSta - 0.5,
+ kNbinsNofSta, -0.5, kNbinsNofSta - 0.5);
+
+ fph_ghost_p->SetTitle("Momentum of ghost track;(p/q)_{reco} [GeV/ec]");
+ fph_ghost_pt->SetTitle("Transverse momentum of ghost track;(p_{T}/q)_{reco} [GeV/ec]");
+ fph_ghost_phi->SetTitle("Azimuthal angle of ghost track;#phi_{reco} [rad]");
+ fph_ghost_nhits->SetTitle("Number of hits in a ghost track;N_{hits}");
+ fph_ghost_fsta->SetTitle("First station of ghost track;ID_{f.st.}");
+ fph_ghost_purity->SetTitle("Purity of ghost track;purity [%]");
+ fph_ghost_chi2_ndf->SetTitle("Fit #chi^{2}/NDF of ghost track;(#chi^{2}/NDF)_{reco}");
+ fph_ghost_prob->SetTitle("Fit probability of ghost track;Prob_{reco}");
+ fph_ghost_tx->SetTitle("Slope along x-axis of ghost tracks;t_{x}");
+ fph_ghost_ty->SetTitle("Slope along y-axis of ghost tracks;t_{y}");
+ fph_ghost_fhitR->SetTitle("Distance of the first hit from z-axis;R [cm]");
+ fph_ghost_nhits_vs_p->SetTitle("Number of hits vs. momentum for ghost track;(p/q)_{reco} [GeV/ec];N_{hits}");
+ fph_ghost_fsta_vs_p->SetTitle("First station vs. momentum for ghost track;(p/q)_{reco} [GeV/ec];ID_{f.st.}");
+ fph_ghost_lsta_vs_fsta->SetTitle("Last station vs. first station for ghost track;ID_{first st.};ID_{last st.}");
+
+ // ----- Reset event counter
+ fNofEvents = 0;
+
+ LOG(info) << "CA QA: Initialized histograms: (directory \033[1;31m" << gDirectory->GetName() << "\033[0m)";
+ for (auto& entry : fmpHistoList) {
+ LOG(info) << "\t- " << entry.first;
+ }
+
+ gDirectory = currentDir;
+}
// ---------------------------------------------------------------------------------------------------------------------
//
-void Qa::FillHistograms() {}
+void Qa::FillHistograms()
+{
+ assert(fpParameters);
+ assert(fpMCData);
+ assert(fpvRecoTracks);
+ assert(fpvHits);
+
+ // Update event counter
+ fNofEvents++;
+
+
+ // ************************************************************************
+ // ** 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
+
+
+ // *************************************
+ // ** Fill distributions of MC-tracks **
+ // *************************************
+
+ //for (const auto& mcTrk: fpMCData->GetTrackContainer()) {
+ //
+ //}
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void Qa::SaveHistograms()
+{
+ fpOutFile->cd();
+ RecursiveWrite(fpHistoDir);
+ LOG(info) << "CA QA: Histograms have been written to \033[1;31m" << fpOutFile->GetName() << "\033[0m";
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void Qa::RecursiveWrite(TObject* pObj)
+{
+ if (!pObj->IsFolder()) { pObj->Write(); }
+ else {
+ TDirectory* pObjDir = dynamic_cast<TDirectory*>(pObj);
+ LOG_IF(fatal, !pObjDir) << "CA QA: Incorrect dynamic cast from the input pointer to RecursiveWrite function";
+
+ TDirectory* pCurDir = gDirectory;
+ TDirectory* pSubDir = pCurDir->mkdir(pObjDir->GetName());
+ pSubDir->cd();
+ TList* pObjList = pObjDir->GetList();
+ TIter it(pObjList);
+ while (TObject* pNext = it()) {
+ RecursiveWrite(pNext);
+ }
+ pCurDir->cd();
+ }
+}
// ---------------------------------------------------------------------------------------------------------------------
//
-void Qa::SaveHistograms() {}
+void Qa::SetOutputFilename(const char* filename)
+{
+ fpOutFile = std::make_unique<TFile>(filename, "RECREATE");
+ if (!fpOutFile.get() || fpOutFile->IsZombie()) {
+ LOG(fatal) << "CA QA: ROOT histogram output file \"" << filename << "\" cannot be created";
+ }
+ fpHistoDir = gROOT->mkdir("histograms");
+}
diff --git a/reco/L1/catools/CaToolsQa.h b/reco/L1/catools/CaToolsQa.h
index 87a58b31b9..e6c05fed33 100644
--- a/reco/L1/catools/CaToolsQa.h
+++ b/reco/L1/catools/CaToolsQa.h
@@ -1,6 +1,6 @@
-/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+/* Copyright (C) 2022-2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
SPDX-License-Identifier: GPL-3.0-only
- Authors: Sergey Gorbunov, Sergei Zharko [committer] */
+ Authors: Sergei Zharko [committer] */
/// \file CaToolsQa.h
/// \brief Tracking performance class (header)
@@ -10,14 +10,22 @@
#ifndef CaToolsQa_h
#define CaToolsQa_h 1
+#include "Logger.h"
+
+#include "TFile.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TProfile.h"
+#include "TROOT.h"
+
+#include <string_view>
+
#include "L1Vector.h"
class L1Parameters;
class TDirectory;
-class TH1;
-class TH1F;
-class TH2F;
-class TProfile;
+class CbmL1HitDebugInfo;
+class CbmL1Track;
namespace ca::tools
{
@@ -54,7 +62,7 @@ namespace ca::tools
void InitHistograms();
/// Fills histograms
- /// \note should be called in the end of the run
+ /// \note should be called in the end of the event
void FillHistograms();
/// Saves histograms to the file
@@ -65,28 +73,176 @@ namespace ca::tools
// ** Accessors **
// ***********************
+ /// Sets pointer to hit container from reference
+ void SetHitContainer(const L1Vector<CbmL1HitDebugInfo>& vHits) { fpvHits = &vHits; }
+
+ /// Sets pointer to hit container from pointer
+ void SetHitContainer(const L1Vector<CbmL1HitDebugInfo>* pvHits)
+ {
+ assert(pvHits);
+ fpvHits = pvHits;
+ }
+
+
/// Sets pointer to MC data object
void SetMCData(const MCData* pMCData) { fpMCData = pMCData; }
+ /// Sets output filename
+ void SetOutputFilename(const char* filename);
+
/// Sets pointer to tracking parameters object
void SetParameters(const L1Parameters* pParameters) { fpParameters = pParameters; }
+ /// Sets pointer to the vector of reconstructed tracks from reference
+ void SetRecoTrackContainer(const L1Vector<CbmL1Track>& vRecoTracks) { fpvRecoTracks = &vRecoTracks; }
+
+ /// Sets pointer to the vector of reconstructed tracks from pointer
+ void SetRecoTrackContainer(const L1Vector<CbmL1Track>* pvRecoTracks)
+ {
+ assert(pvRecoTracks);
+ fpvRecoTracks = pvRecoTracks;
+ }
+
private:
- const L1Parameters* fpParameters = nullptr; ///< Instance of the tracking core class parameters
- const MCData* fpMCData = nullptr; ///< Container for MC information on the event
+ /// Function, which creates a histogram and initializes it
+ /// \tparam T Type of the histogram (\note T should not be a pointer)
+ /// \param name Name of histogram
+ /// \param args The rest of the arguments, which will be passed to the histogram constructor
+ template<typename T, typename... Args>
+ T* MakeHisto(const char* name, Args... args);
+
+ /// Writes recursively all registered objects inside a current directory
+ /// \param pObj Pointer to object to write
+ void RecursiveWrite(TObject* pObj);
+
+ const L1Parameters* fpParameters = nullptr; ///< Instance of the tracking core class parameters
+ const MCData* fpMCData = nullptr; ///< Container for MC information on the event
+ const L1Vector<CbmL1Track>* fpvRecoTracks = nullptr; ///< Pointer to external reconstructed track container
+ const L1Vector<CbmL1HitDebugInfo>* fpvHits = nullptr; ///< Pointer to external hit container
// ***********************
// ** Histograms **
// ***********************
- TDirectory* fpHistoDirectory = nullptr; ///< directory object to store performance histograms
-
- TH1F* fphRegMcMom = nullptr;
- TH1F* fphAccMCMom = nullptr;
-
- L1Vector<TH1*> fvpHistoList = {"ca::tools::Qa::fvpHistoList"}; ///< List of histograms for iteration over them
+ TDirectory* fpHistoDir = nullptr; ///< directory to store QA histograms
+
+ std::unique_ptr<TFile> fpOutFile = nullptr;
+
+ // ** Number of bins for different quantities
+
+ static constexpr int kNbinsRecoP = 50; ///< momentum of reconstructed tracks
+ static constexpr int kNbinsRecoPt = 50; ///< transverse momentum of reconstructed tracks
+ static constexpr int kNbinsRecoPhi = 68; ///< azimuthal angle of reconstructed tracks
+ static constexpr int kNbinsNofHits = 10; ///< number of hits of reconstructed tracks
+ static constexpr int kNbinsPurity = 100; ///< track purity
+ static constexpr int kNbinsChi2NDF = 50; ///< fit chi-square over NDF
+ static constexpr int kNbinsProb = 505; ///< fit probability
+ static constexpr int kNbinsNofSta = 15; ///< number of stations
+ static constexpr int kNbinsTx = 50; ///< slope along x-axis
+ static constexpr int kNbinsTy = 50; ///< slope along y-axis
+ static constexpr int kNbinsHitDist = 50; ///< distance from transverse hit position to z-axis
+
+ // ** Upper bound for quantity range for histograms
+ static constexpr double kAxMaxRecoP = 5.; ///< momentum of reconstructed tracks [GeV/c]
+ static constexpr double kAxMaxRecoPt = 5.; ///< transverse momentum of reconstructed tracks [GeV/c]
+ static constexpr double kAxMaxRecoPhi = 3.2; ///< azimuthal angle of reconstructed tracks [rad]
+ static constexpr double kAxMaxChi2NDF = 100.; ///< fit chi-square over NDF
+
+
+ // ** 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_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
+
+ // TODO: Add other histograms
+
+ // Registered primary tracks
+ // Registered secondary tracks
+ // Reconstructable primary tracks
+ // Reconstructable secondary tracks
+ // Reconstructed primary tracks
+ // Reconstructed secondary tracks
+ // Reconstructed tracks
+ // Ghost tracks
+
+ // ** Histograms appearance attributes
+ static constexpr double kHtextSize = 0.04;
+
+
+ /// A map of the histogram pointers. The key is a name of the histogram
+ std::unordered_map<std::string_view, TH1*> fmpHistoList; ///< List of pointers to histograms
+
+ // *******************************
+ // ** Utility variables **
+ // *******************************
+
+ int fNofEvents = 0; ///< Event counter
};
} // namespace ca::tools
+// **********************************************************
+// ** Inline and template function implementations **
+// **********************************************************
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+template<typename T, typename... Args>
+T* ca::tools::Qa::MakeHisto(const char* name, Args... args)
+{
+ // Check, if the histogram with a given name was already created. If so, delete it
+ if (gROOT->FindObject(name)) {
+ LOG(warn) << "CA QA: A histogram with name \"" << name << "\" was previously created and will be deleted now "
+ << "to avoid memory leaks";
+ T* pHisto = (T*) gROOT->FindObject(name);
+ delete pHisto;
+ }
+
+ // Create a new histogram or profile
+ T* pHisto = new T(name, args...);
+
+ // Register histogram in the current TDirectory
+ pHisto->SetDirectory(fpHistoDir);
+
+ // Register histogram in the list
+ fmpHistoList[pHisto->GetName()] = (TH1*) pHisto;
+
+ pHisto->Sumw2();
+
+ // Appearance settings
+ std::array<TAxis*, 2> vpAxes = {pHisto->GetXaxis(), pHisto->GetYaxis()};
+ for (auto* pAxis : vpAxes) {
+ pAxis->SetTitleSize(kHtextSize);
+ pAxis->SetLabelSize(kHtextSize);
+ }
+
+ return pHisto;
+}
+
+
#endif // CaToolsQa_h
--
GitLab