From 10d7b9d11e739e7f59c0ad317183bf9adbb51c97 Mon Sep 17 00:00:00 2001
From: "s.zharko@gsi.de" <s.zharko@gsi.de>
Date: Fri, 3 Feb 2023 19:24:46 +0100
Subject: [PATCH] CA: QA task for tracking input from MuCh, TRD and TOF (v.2)
CbmTofInteraction class is introduced to address the problem of multiple MC points for one interaction.
---
core/detectors/tof/CbmTofTrackingInterface.h | 2 +-
core/qa/CMakeLists.txt | 3 +-
core/qa/CbmQaEff.cxx | 63 +
core/qa/CbmQaEff.h | 73 ++
core/qa/CbmQaTable.cxx | 63 +-
core/qa/CbmQaTable.h | 30 +-
core/qa/CbmQaTask.cxx | 24 +-
core/qa/CbmQaTask.h | 98 +-
macro/run/run_reco_qa.C | 19 +-
reco/L1/CMakeLists.txt | 8 +-
reco/L1/CbmL1.cxx | 3 +-
reco/L1/CbmL1ReadEvent.cxx | 2 +-
reco/L1/L1Algo/L1Constants.h | 7 +
reco/L1/L1Algo/L1Undef.h | 70 ++
reco/L1/L1LinkDef.h | 2 +
reco/L1/catools/CaToolsLinkKey.h | 52 +
reco/L1/qa/CbmCaInputQaBase.h | 80 ++
reco/L1/qa/CbmCaInputQaMuch.cxx | 1020 ++++++++++++++++
reco/L1/qa/CbmCaInputQaMuch.h | 166 +++
reco/L1/qa/CbmCaInputQaSts.cxx | 200 ++--
reco/L1/qa/CbmCaInputQaSts.h | 41 +-
reco/L1/qa/CbmCaInputQaTof.cxx | 1108 ++++++++++++++++++
reco/L1/qa/CbmCaInputQaTof.h | 178 +++
reco/L1/qa/CbmCaInputQaTrd.cxx | 1056 +++++++++++++++++
reco/L1/qa/CbmCaInputQaTrd.h | 167 +++
reco/L1/qa/CbmCaQaCuts.h | 3 +-
reco/L1/qa/CbmTofInteraction.cxx | 99 ++
reco/L1/qa/CbmTofInteraction.h | 77 ++
28 files changed, 4509 insertions(+), 205 deletions(-)
create mode 100644 core/qa/CbmQaEff.cxx
create mode 100644 core/qa/CbmQaEff.h
create mode 100644 reco/L1/L1Algo/L1Undef.h
create mode 100644 reco/L1/catools/CaToolsLinkKey.h
create mode 100644 reco/L1/qa/CbmCaInputQaBase.h
create mode 100644 reco/L1/qa/CbmCaInputQaMuch.cxx
create mode 100644 reco/L1/qa/CbmCaInputQaMuch.h
create mode 100644 reco/L1/qa/CbmCaInputQaTof.cxx
create mode 100644 reco/L1/qa/CbmCaInputQaTof.h
create mode 100644 reco/L1/qa/CbmCaInputQaTrd.cxx
create mode 100644 reco/L1/qa/CbmCaInputQaTrd.h
create mode 100644 reco/L1/qa/CbmTofInteraction.cxx
create mode 100644 reco/L1/qa/CbmTofInteraction.h
diff --git a/core/detectors/tof/CbmTofTrackingInterface.h b/core/detectors/tof/CbmTofTrackingInterface.h
index 205e29599b..47b0b37edd 100644
--- a/core/detectors/tof/CbmTofTrackingInterface.h
+++ b/core/detectors/tof/CbmTofTrackingInterface.h
@@ -108,7 +108,7 @@ public:
// CbmTofAddress::GetSmId(address),
// CbmTofAddress::GetRpcId(address));
// NOTE: Implement, when the "mcbm_beam_2021_07_surveyed" parameters will be fixed
-
+ // TODO: Invesitgate problem in mcbm_beam_2021_07_surveyed
return -1; // iSt;
}
diff --git a/core/qa/CMakeLists.txt b/core/qa/CMakeLists.txt
index 3589694a85..f59b755abe 100644
--- a/core/qa/CMakeLists.txt
+++ b/core/qa/CMakeLists.txt
@@ -4,6 +4,7 @@ set(INCLUDE_DIRECTORIES
set(SRCS
CbmQaCanvas.cxx
+ CbmQaEff.cxx
CbmQaPie.cxx
CbmQaHist.cxx
CbmQaTable.cxx
@@ -27,6 +28,6 @@ set(INTERFACE_DEPENDENCIES
generate_cbm_library()
-Install(FILES CbmQaTask.h CbmQaCanvas.h CbmQaTable.h CbmQaHist.h
+Install(FILES CbmQaTask.h CbmQaCanvas.h CbmQaTable.h CbmQaHist.h CbmQaEff.h CbmQaPie.h
DESTINATION include
)
diff --git a/core/qa/CbmQaEff.cxx b/core/qa/CbmQaEff.cxx
new file mode 100644
index 0000000000..02d295c68c
--- /dev/null
+++ b/core/qa/CbmQaEff.cxx
@@ -0,0 +1,63 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// \file CbmQaEff.cxx
+/// \date 20.01.2023
+/// \author S.Zharko <s.zharko@gsi.de>
+/// \brief Implementation of CbmQaEff class
+
+#include "CbmQaEff.h"
+
+#include "Logger.h"
+
+#include "TString.h"
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+CbmQaEff::CbmQaEff() : TEfficiency() {}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+CbmQaEff::CbmQaEff(const CbmQaEff& other) : TEfficiency(other) {}
+
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+CbmQaEff* CbmQaEff::Integral(double lo, double hi)
+{
+ // Underlying histograms with passed and total events
+ auto* pPassed = (TH1D*) (this->GetPassedHistogram());
+ auto* pTotal = (TH1D*) (this->GetTotalHistogram());
+
+ if (!pPassed) { return nullptr; }
+
+ // Integration range
+ double range[2] = {0};
+ if (lo < hi) {
+ range[0] = lo;
+ range[1] = hi;
+ }
+ else {
+ range[0] = pPassed->GetXaxis()->GetXmin();
+ range[1] = pPassed->GetXaxis()->GetXmax();
+ }
+
+ // Re-binned histograms for a selected integration range
+ auto& histPassedReb = *(pPassed->Rebin(1, "ptmp", range));
+ auto& histTotalReb = *(pTotal->Rebin(1, "ttmp", range));
+
+ LOG(info) << "DEBUG: " << this->GetName() << ": passed: " << pPassed->GetEntries();
+ LOG(info) << "DEBUG: " << this->GetName() << ": total: " << pTotal->GetEntries();
+ LOG(info) << "DEBUG: " << this->GetName() << ": passed: " << histPassedReb.GetBinContent(1) << ' '
+ << histPassedReb.GetBinError(1);
+ LOG(info) << "DEBUG: " << this->GetName() << ": total: " << histTotalReb.GetBinContent(1) << ' '
+ << histPassedReb.GetBinError(1);
+
+ TString sName = Form("%s_integrated", this->GetName());
+
+ // New efficiency
+ auto* pIntEff = new CbmQaEff(TEfficiency(histPassedReb, histTotalReb));
+ pIntEff->SetName(sName);
+ return pIntEff;
+}
diff --git a/core/qa/CbmQaEff.h b/core/qa/CbmQaEff.h
new file mode 100644
index 0000000000..40977f08a3
--- /dev/null
+++ b/core/qa/CbmQaEff.h
@@ -0,0 +1,73 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// \file CbmQaEff.h
+/// \date 20.01.2023
+/// \author S.Zharko <s.zharko@gsi.de>
+/// \brief Declaration of CbmQaEff class
+
+
+#ifndef CbmQaEff_h
+#define CbmQaEff_h 1
+
+#include "CbmQaCanvas.h"
+
+#include "Logger.h"
+
+#include "TEfficiency.h"
+#include "TFitResultPtr.h"
+#include "TGraphAsymmErrors.h"
+#include "TH2.h"
+#include "TPaveStats.h"
+#include "TStyle.h"
+#include "TVirtualPad.h"
+
+
+/// Implementation of ROOT TEfficiency class, which adds handy functionality and improves fitting and drawing
+///
+class CbmQaEff : public TEfficiency {
+public:
+ /// Default constructor
+ CbmQaEff();
+
+ /// Copy constructor
+ CbmQaEff(const CbmQaEff& other);
+
+ /// Other constructors
+ /// \tparam Args Variadic template parameter, which is passed to the constructor of base class
+ template<typename... Args>
+ CbmQaEff(Args... args);
+
+ /// Destructor
+ ~CbmQaEff() = default;
+
+ /// Fit method reimplementation
+ template<typename... Args>
+ TFitResultPtr Fit(Args... args)
+ {
+ return TEfficiency::Fit(args...);
+ }
+
+ /// Integrates efficiency over the range
+ /// \param lo Lower bound of integration range
+ /// \param hi Higher bound of integration range
+ /// \return Pointer to efficiency object, which contains only one point
+ CbmQaEff* Integral(double lo, double hi);
+
+
+ ClassDef(CbmQaEff, 1);
+};
+
+// **********************************************************
+// ** Template and inline functions implementation **
+// **********************************************************
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+template<typename... Args>
+CbmQaEff::CbmQaEff(Args... args) : TEfficiency(args...)
+{
+}
+
+#endif // CbmQaEff_h
diff --git a/core/qa/CbmQaTable.cxx b/core/qa/CbmQaTable.cxx
index 9e284febdc..d7f4859918 100644
--- a/core/qa/CbmQaTable.cxx
+++ b/core/qa/CbmQaTable.cxx
@@ -31,7 +31,7 @@ CbmQaTable::CbmQaTable(const char* name, const char* title, Int_t nRows, Int_t n
//
// Setup default style of the table
TH2D::SetStats(kFALSE);
- TH2D::SetOption("text");
+ TH2D::SetOption("textX+");
TH2D::GetXaxis()->SetTickLength(0.);
TH2D::GetYaxis()->SetTickLength(0.);
TH2D::GetXaxis()->SetLabelFont(kDefaultFontStyle);
@@ -39,7 +39,8 @@ CbmQaTable::CbmQaTable(const char* name, const char* title, Int_t nRows, Int_t n
// Define basic names of rows and columns
for (Int_t iRow = 1; iRow <= fNrows; ++iRow) {
- TH2D::GetYaxis()->SetBinLabel(fNrows - iRow + 1, TString::Format("row %d", iRow).Data());
+ //TH2D::GetYaxis()->SetBinLabel(fNrows - iRow + 1, TString::Format("row %d", iRow).Data());
+ TH2D::GetYaxis()->SetBinLabel(fNrows - iRow + 1, "");
}
for (Int_t iCol = 1; iCol <= fNcols; ++iCol) {
TH2D::GetXaxis()->SetBinLabel(iCol, TString::Format("col %d", iCol).Data());
@@ -90,10 +91,35 @@ void CbmQaTable::SetNamesOfRows(const std::vector<std::string>& names)
//
//------------------------------------------------------------------------------------------------------------------------
//
-std::string CbmQaTable::ToString() const
+std::string CbmQaTable::ToString(int prec) const
{
std::stringstream aStream;
- aStream << (*this);
+ aStream.setf(std::ios::fixed);
+ aStream.setf(std::ios::showpoint);
+ aStream.setf(std::ios::left);
+ aStream.precision(prec);
+
+ aStream << "Table: " << GetTitle() << '\n';
+ aStream << std::setw(fColWidth) << std::setfill(' ') << ' ' << ' ';
+
+ for (Int_t iCol = 0; iCol < fNcols; ++iCol) {
+ std::string entry = std::string(this->GetXaxis()->GetBinLabel(iCol + 1));
+ if (Int_t(entry.size()) + 3 > fColWidth) { entry = entry.substr(0, fColWidth - 3) + "..."; }
+ aStream << std::setw(fColWidth) << std::setfill(' ') << entry << ' ';
+ }
+ aStream << '\n';
+
+ for (Int_t iRow = 0; iRow < fNrows; ++iRow) {
+ // Print row title
+ std::string entry = std::string(GetYaxis()->GetBinLabel(fNrows - iRow));
+ if (static_cast<Int_t>(entry.size()) > fColWidth) { entry = entry.substr(0, fColWidth - 3) + "..."; }
+ aStream << std::setw(fColWidth) << std::setfill(' ') << entry << ' ';
+
+ for (Int_t iCol = 0; iCol < fNcols; ++iCol) {
+ aStream << std::setw(fColWidth) << std::setfill(' ') << GetCell(iRow, iCol) << ' ';
+ }
+ aStream << '\n';
+ }
return aStream.str();
}
//
@@ -120,34 +146,7 @@ void CbmQaTable::SetTextSize(Float_t size)
//
std::ostream& operator<<(std::ostream& out, const CbmQaTable& aTable)
{
- out.setf(std::ios::fixed);
- out.setf(std::ios::showpoint);
- out.precision(3);
- out.setf(std::ios::left);
- // Print column titles
- out << std::setw(CbmQaTable::kRowTitlesSetwPar) << std::setfill(' ') << ' ' << ' '; // top-left cell, always
- for (Int_t iCol = 1; iCol <= aTable.fNcols; ++iCol) {
- std::string entry = std::string(aTable.GetXaxis()->GetBinLabel(iCol));
- if (static_cast<Int_t>(entry.size()) > CbmQaTable::kDefaultSetwPar) {
- entry = entry.substr(0, CbmQaTable::kDefaultSetwPar - 3) + "...";
- }
- out << std::setw(CbmQaTable::kDefaultSetwPar) << std::setfill(' ') << entry << ' ';
- }
- out << '\n';
-
- for (Int_t iRow = 0; iRow < aTable.fNrows; ++iRow) {
- // Print row title
- std::string entry = std::string(aTable.GetYaxis()->GetBinLabel(aTable.fNrows - iRow));
- if (static_cast<Int_t>(entry.size()) > CbmQaTable::kRowTitlesSetwPar) {
- entry = entry.substr(0, CbmQaTable::kDefaultSetwPar - 3) + "...";
- }
- out << std::setw(CbmQaTable::kRowTitlesSetwPar) << std::setfill(' ') << entry << ' ';
-
- for (Int_t iCol = 0; iCol < aTable.fNcols; ++iCol) {
- out << std::setw(CbmQaTable::kDefaultSetwPar) << std::setfill(' ') << aTable.GetCell(iRow, iCol) << ' ';
- }
- out << '\n';
- }
+ out << aTable.ToString(3);
return out;
}
//
diff --git a/core/qa/CbmQaTable.h b/core/qa/CbmQaTable.h
index 97115ff8db..caba486994 100644
--- a/core/qa/CbmQaTable.h
+++ b/core/qa/CbmQaTable.h
@@ -18,47 +18,57 @@
#include <fstream>
#include <string>
-// TODO: We need a method, which sets and gets cell values! (S.Zharko)
+
+/// TODO: SZh, 30.01.2023: Override THistPainter::PaintText() to add zeroes in tables
class CbmQaTable : public TH2D {
public:
- //
- // CONSTRUCTORS AND DESTRUCTORS
- //
/// Default constructor
CbmQaTable() : TH2D() {}
+
/// Constructor from number of rows and columns
CbmQaTable(const char* name, const char* title, Int_t nRows, Int_t nCols);
+
/// Destructor
virtual ~CbmQaTable();
/// Dumps table content into a string
- std::string ToString() const;
+ /// \param prec Precision of numbers
+ std::string ToString(int prec) const;
+
/// Dumps table content into a text file. File open mode is also controllable, for example, use
/// mode = std::ios_base::app to append the table into an existing file
void ToTextFile(const std::string& fileName, std::ios_base::openmode mode = std::ios_base::out) const;
- //
- // GETTERS
- //
/// Gets cell content. Please mind, that the signature and result of this function is different to TH2D::GetBinContent
Double_t GetCell(Int_t iRow, Int_t iCol) const;
+
/// Gets cell error. Please mind, that the signature and result of this function is different to TH2D::GetBinError
Double_t GetCellError(Int_t iRow, Int_t iCol) const;
+
/// Sets number of rows
Int_t GetNrows() const { return fNrows; }
+
/// Sets number of columns
Int_t GetNcols() const { return fNcols; }
+
/// Sets cell content and error. Please mind, that the signature and result of this function
/// is different to TH2D::SetBinContent and TH2D::SetBinError
void SetCell(Int_t iRow, Int_t iCol, Double_t content, Double_t error = 0.);
+
/// Sets the names of table columns
void SetNamesOfCols(const std::vector<std::string>& names);
+
/// Sets the names of table rows
void SetNamesOfRows(const std::vector<std::string>& names);
+
/// Sets size of the text
void SetTextSize(Float_t size = 0.03);
+ /// Sets width of column in log output
+ /// \param width Width of column [number of symbols]
+ void SetColWidth(Int_t width) { fColWidth = width; }
+
/// Dumps table content into a stream
friend std::ostream& operator<<(std::ostream& out, const CbmQaTable& aTable);
@@ -75,11 +85,13 @@ private:
Int_t fNcols {0}; ///< number of columns in a table
Int_t fNrows {0}; ///< number of rows in a table
+ Int_t fColWidth {12}; ///< Width of column in number of symbols
+
// Some hard-coded constants
static constexpr Float_t kDefaultTextSize {0.03}; ///< default size of text
static constexpr Style_t kDefaultFontStyle {62}; ///< default text style
- static constexpr Int_t kDefaultSetwPar {12}; ///< default size of entry in std::setw() for columns
+ static constexpr Int_t kDefaultSetwPar {20}; ///< default size of entry in std::setw() for columns
static constexpr Int_t kRowTitlesSetwPar {30}; ///< size of entry in std::setw() for row titles
static constexpr Int_t kValuesPrecision {3}; ///< precision of output parameters
// TODO: Apply this precision and other options to graphical verion of the table
diff --git a/core/qa/CbmQaTask.cxx b/core/qa/CbmQaTask.cxx
index d74cd4582c..78ce9c0386 100644
--- a/core/qa/CbmQaTask.cxx
+++ b/core/qa/CbmQaTask.cxx
@@ -1,4 +1,4 @@
-/* opyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
SPDX-License-Identifier: GPL-3.0-only
Authors: Sergei Zharko [committer] */
@@ -24,7 +24,12 @@ ClassImp(CbmQaTask);
// ---------------------------------------------------------------------------------------------------------------------
//
-CbmQaTask::CbmQaTask(const char* name, int verbose, bool isMCUsed) : FairTask(name, verbose), fbUseMC(isMCUsed) {}
+CbmQaTask::CbmQaTask(const char* name, const char* prefix, int verbose, bool isMCUsed)
+ : FairTask(name, verbose)
+ , fsPrefix(prefix)
+ , fbUseMC(isMCUsed)
+{
+}
// ---------------------------------------------------------------------------------------------------------------------
//
@@ -40,10 +45,9 @@ void CbmQaTask::Exec(Option_t* /*option*/)
void CbmQaTask::Finish()
{
// Processes histograms in the end of the run
- AnalyzeHistograms();
-
- // Check QA
- Check();
+ LOG(info) << fName << ": Checking QA...";
+ bool qaResult = Check();
+ LOG(info) << fName << ": Checking QA: " << (qaResult ? "\033[1;32mpassed\033[0m" : "\033[1;31mfailed\033[0m");
// Processes canvases in the end of the run
LOG_IF(info, fVerbose > 1) << fName << ": initializing canvases";
@@ -52,7 +56,6 @@ void CbmQaTask::Finish()
// Write the root folder to sinker
FairSink* pSink = FairRootManager::Instance()->GetSink();
LOG_IF(fatal, !pSink) << fName << ": sink file was not found";
-
pSink->WriteObject(fpFolderRoot.get(), nullptr);
}
@@ -150,13 +153,6 @@ void CbmQaTask::FillHistograms()
LOG_IF(info, fVerbose > 1) << fName << ": histogram filling function is not defined";
}
-// ---------------------------------------------------------------------------------------------------------------------
-//
-void CbmQaTask::AnalyzeHistograms()
-{
- LOG_IF(info, fVerbose > 1) << fName << ": no action on histograms in the end of the run is provided";
-}
-
// ---------------------------------------------------------------------------------------------------------------------
//
InitStatus CbmQaTask::InitCanvases()
diff --git a/core/qa/CbmQaTask.h b/core/qa/CbmQaTask.h
index c18e596485..2415e28a41 100644
--- a/core/qa/CbmQaTask.h
+++ b/core/qa/CbmQaTask.h
@@ -10,6 +10,8 @@
#ifndef CbmQaTask_h
#define CbmQaTask_h 1
+#include "CbmQaTable.h"
+
#include "FairTask.h"
#include "Logger.h"
@@ -20,6 +22,7 @@
#include "TParameter.h"
#include "TROOT.h"
+#include <algorithm>
#include <string_view>
/// Class CbmQaTask is to be inherited with a particular QA-task. It provides mechanisms for storage and management
@@ -28,9 +31,10 @@ class CbmQaTask : public FairTask {
public:
/// Constructor from parameters
/// \param name Name of the task
+ /// \param prefix Unique prefix for writeable root objects
/// \param verbose Verbose level
/// \param isMCUsed Flag: true - MC information is used, false - only reconstructed data QA is processed
- CbmQaTask(const char* name, int verbose, bool isMCUsed);
+ CbmQaTask(const char* name, const char* prefix, int verbose, bool isMCUsed);
/// Default constructor
CbmQaTask() = delete; // TODO: Let's see, what can happen, if one deletes default constructor
@@ -66,9 +70,8 @@ protected:
// ** Functions accessible inside the derived classes **
// *****************************************************
- /// Method to check, if the QA results are acceptable
- /// NOTE: Here one can add a geometry check ...
- virtual bool Check() const { return true; }
+ /// \brief Method to check, if the QA results are acceptable
+ virtual bool Check() = 0;
/// De-initialize the task
virtual void DeInit();
@@ -85,30 +88,33 @@ protected:
/// Method to fill histograms per event or time-slice
virtual void FillHistograms();
- /// Method to process histograms in the end of the run (optional)
- virtual void AnalyzeHistograms();
-
- /// Creates, initializes and registers the histogram
+ /// Creates, initializes and registers a canvas
/// \tparam T Type of the canvas: TCanvas or CbmQaCanvas (\note T should not be a pointer)
/// \param name Name of canvas
/// \param args The rest of the arguments, which will be passed to the canvas constructor
template<typename T, typename... Args>
T* MakeCanvas(const char* name, Args... args);
- /// Creates, initializes and registers the efficiency
+ /// Creates, initializes and registers an efficiency
/// \tparam T Type of the canvas: either TProfile, TProfile2D or TEfficiency (\note T should not be a pointer)
/// \param name Name of canvas
/// \param args The rest of the arguments, which will be passed to the efficiency constructor
template<typename T, typename... Args>
T* MakeEfficiency(const char* name, Args... args);
- /// Creates, initializes and registers the histogram
+ /// Creates, initializes and registers a histogram
/// \tparam T Type of the histogram (\note T should not be a pointer)
/// \param name Name of histogram
/// \param args The rest of the arguments, which will be passed to the histogram constructor
template<typename T, typename... Args>
T* MakeHisto(const char* name, Args... args);
+ /// 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);
@@ -123,6 +129,20 @@ protected:
int GetEventNumber() const { return fNofEvents.GetVal(); }
+ // ***********************
+ // ** Utility functions **
+ // ***********************
+
+ /// Checks range of variable
+ /// \param name Name of the variable
+ /// \param var Variable to check
+ /// \param lo Lower limit of the variable
+ /// \param hi Upper limit of the variable
+ /// \return False, if variable exits the range
+ template<typename T>
+ bool CheckRange(const std::string_view& name, const T& var, const T& lo, const T& hi) const;
+
+
private:
// **********************************************************
// ** Functions accessible inside the CbmQaTask class only **
@@ -135,14 +155,15 @@ private:
// ** Private variables **
// ***********************
- // Flags
- bool fbUseMC = false; ///< Flag, if MC is used
-
// I/O
std::unique_ptr<TFolder> fpFolderRoot = nullptr; ///< Root folder to store histograms and canvases
TFolder* fpFolderHist = nullptr; ///< Folder for raw histograms
TFolder* fpFolderCanv = nullptr; ///< Folder for canvases
TFolder* fpFolderEff = nullptr; ///< Folder for efficiencies
+ TFolder* fpFolderTable = nullptr; ///< Folder for tables
+
+ TString fsPrefix = ""; ///< Unique prefix for all writeable root objects to avoid collisions between different tasks
+ bool fbUseMC = false; ///< Flag, if MC is used
TParameter<int> fNofEvents {"nEvents", 0}; ///< Number of processed events
};
@@ -152,11 +173,28 @@ private:
// ** Inline and template function implementation **
// *************************************************
+// ---------------------------------------------------------------------------------------------------------------------
+//
+template<typename T>
+bool CbmQaTask::CheckRange(const std::string_view& name, const T& var, const T& lo, const T& hi) const
+{
+ if (std::clamp(var, lo, hi) == lo) {
+ LOG(error) << fName << ": " << name << " is found to be under the lower limit (" << var << " < " << lo << ')';
+ return false;
+ }
+ if (std::clamp(var, lo, hi) == hi) {
+ LOG(error) << fName << ": " << name << " is found to be above the upper limit (" << var << " > " << hi << ')';
+ return false;
+ }
+ return true;
+}
+
// ---------------------------------------------------------------------------------------------------------------------
//
template<typename T, typename... Args>
-T* CbmQaTask::MakeCanvas(const char* name, Args... args)
+T* CbmQaTask::MakeCanvas(const char* nameBase, Args... args)
{
+ TString name = fsPrefix + "_" + nameBase;
if (gROOT->FindObject(name)) {
LOG(warn) << fName << ": A canvas with name \"" << name << "\" was previously created and will be deleted now "
<< "to avoid memory leaks";
@@ -180,8 +218,9 @@ T* CbmQaTask::MakeCanvas(const char* name, Args... args)
// ---------------------------------------------------------------------------------------------------------------------
//
template<typename T, typename... Args>
-T* CbmQaTask::MakeEfficiency(const char* name, Args... args)
+T* CbmQaTask::MakeEfficiency(const char* nameBase, Args... args)
{
+ TString name = fsPrefix + "_" + nameBase;
if (gROOT->FindObject(name)) {
LOG(warn) << fName << ": An efficiency with name \"" << name << "\" was previously created and will be deleted now "
<< "to avoid memory leaks";
@@ -203,8 +242,9 @@ T* CbmQaTask::MakeEfficiency(const char* name, Args... args)
// ---------------------------------------------------------------------------------------------------------------------
//
template<typename T, typename... Args>
-T* CbmQaTask::MakeHisto(const char* name, Args... args)
+T* CbmQaTask::MakeHisto(const char* nameBase, Args... args)
{
+ TString name = fsPrefix + "_" + nameBase;
// Check, if the histogram with a given name was already created. If so, delete it
if (gROOT->FindObject(name)) {
LOG(warn) << fName << ": A histogram with name \"" << name << "\" was previously created and will be deleted now "
@@ -215,7 +255,6 @@ T* CbmQaTask::MakeHisto(const char* name, Args... args)
// Create a new histogram or profile
T* pHist = new T(name, args...);
- pHist->Sumw2();
// Register histogram in the folder
if (!fpFolderHist) { fpFolderHist = fpFolderRoot->AddFolder("histograms", "Histograms"); }
@@ -226,4 +265,29 @@ T* CbmQaTask::MakeHisto(const char* name, Args... args)
return pHist;
}
+// ---------------------------------------------------------------------------------------------------------------------
+//
+template<typename... Args>
+CbmQaTable* CbmQaTask::MakeTable(const char* nameBase, Args... args)
+{
+ TString name = fsPrefix + "_" + nameBase;
+ // Check, if the table with a given name was already created. If so, delete it
+ if (gROOT->FindObject(name)) {
+ LOG(warn) << fName << ": A table with name \"" << name << "\" was previously created and will be deleted now "
+ << "to avoid memory leaks";
+ CbmQaTable* pTable = (CbmQaTable*) gROOT->FindObject(name);
+ delete pTable;
+ }
+
+ // Create a new table
+ CbmQaTable* pTable = new CbmQaTable(name, args...);
+
+ // Register table in folder
+ if (!fpFolderTable) { fpFolderTable = fpFolderRoot->AddFolder("tables", "Tables"); }
+ fpFolderTable->Add(pTable);
+
+ return pTable;
+}
+
+
#endif // CbmQaTask_h
diff --git a/macro/run/run_reco_qa.C b/macro/run/run_reco_qa.C
index f286a031ea..446121bace 100644
--- a/macro/run/run_reco_qa.C
+++ b/macro/run/run_reco_qa.C
@@ -187,6 +187,10 @@ void run_reco_qa(TString dataTra = "data/sis100_muon_jpsi_test", TString dataRaw
CbmMuchHitFinderQa* muchHitFinderQa = new CbmMuchHitFinderQa();
muchHitFinderQa->SetGeoFileName(muchParFile);
run->AddTask(muchHitFinderQa);
+
+ auto* pInputQaMuch = new CbmCaInputQaMuch(verbose, bUseMC);
+ pInputQaMuch->SetEfficiencyThrsh(0.5, 0, 100);
+ run->AddTask(new CbmCaInputQaMuch(verbose, bUseMC));
}
// ----- TRD QA ---------------------------------
@@ -198,10 +202,14 @@ void run_reco_qa(TString dataTra = "data/sis100_muon_jpsi_test", TString dataRaw
run->AddTask(new CbmTrdHitProducerQa());
run->AddTask(new CbmTrdCalibTracker());
run->AddTask(new CbmTrackerInputQaTrd()); // Tracker requirements to TRD
+
+ auto* pInputQaTrd = new CbmCaInputQaTrd(verbose, bUseMC);
+ pInputQaTrd->SetEfficiencyThrsh(0.5, 0, 100); // Integration over the whole range
+ run->AddTask(new CbmCaInputQaTrd(verbose, bUseMC)); // Tracker requirements to TRD
}
// ------------------------------------------------------------------------
- // ----- TRD QA ---------------------------------
+ // ----- TOF QA ---------------------------------
if (CbmSetup::Instance()->IsActive(ECbmModuleId::kTof)) {
// TODO: SZh 19.10.2022:
// After the proper TOF digi parameters initialization it is appeared,
@@ -212,6 +220,10 @@ void run_reco_qa(TString dataTra = "data/sis100_muon_jpsi_test", TString dataRaw
// resolved.
//run->AddTask(new CbmTrackerInputQaTof()); // Tracker requirements to TRD
+
+ auto* pInputQaTof = new CbmCaInputQaTof(verbose, bUseMC);
+ pInputQaTof->SetEfficiencyThrsh(0.5, 0, 100);
+ run->AddTask(new CbmCaInputQaTof(verbose, bUseMC));
}
// ------------------------------------------------------------------------
@@ -219,8 +231,9 @@ void run_reco_qa(TString dataTra = "data/sis100_muon_jpsi_test", TString dataRaw
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 CbmCaInputQaSts(verbose, bUseMC));
+
+ auto* pInputQaSts = new CbmCaInputQaSts(verbose, bUseMC);
+ run->AddTask(pInputQaSts);
}
// ------------------------------------------------------------------------
diff --git a/reco/L1/CMakeLists.txt b/reco/L1/CMakeLists.txt
index e9c9c00351..2532b316b5 100644
--- a/reco/L1/CMakeLists.txt
+++ b/reco/L1/CMakeLists.txt
@@ -77,8 +77,11 @@ set(SRCS
qa/CbmTrackerInputQaTrd.cxx
qa/CbmTrackerInputQaTof.cxx
qa/CbmTrackingInputQaSts.cxx
- qa/CbmCaInputQaMuch.cxx
qa/CbmCaInputQaSts.cxx
+ qa/CbmCaInputQaMuch.cxx
+ qa/CbmCaInputQaTrd.cxx
+ qa/CbmCaInputQaTof.cxx
+ qa/CbmTofInteraction.cxx # Tests
)
set(NO_DICT_SRCS
@@ -96,8 +99,11 @@ set(HEADERS
CbmL1Vtx.h
L1Algo/L1Def.h
L1Algo/L1Vector.h
+ L1Algo/L1Undef.h
catools/CaToolsWindowFinder.h
+ catools/CaToolsLinkKey.h
qa/CbmCaQaCuts.h
+
)
diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx
index fa44b00e76..d03df39f8c 100644
--- a/reco/L1/CbmL1.cxx
+++ b/reco/L1/CbmL1.cxx
@@ -926,8 +926,9 @@ void CbmL1::Reconstruct(CbmEvent* event)
else {
fvSortedStsHitsIndexes.clear();
fvSortedStsHitsIndexes.reserve(nStsHits);
- for (int i = 0; i < nStsHits; i++)
+ for (int i = 0; i < nStsHits; i++) {
fvSortedStsHitsIndexes.push_back(i);
+ }
}
// -----------------------------------------------------------------------
diff --git a/reco/L1/CbmL1ReadEvent.cxx b/reco/L1/CbmL1ReadEvent.cxx
index 6640dc4d36..a904f2c4e0 100644
--- a/reco/L1/CbmL1ReadEvent.cxx
+++ b/reco/L1/CbmL1ReadEvent.cxx
@@ -1315,7 +1315,7 @@ void CbmL1::ReadEvent(float& TsStart, float& TsLength, float& /*TsOverlap*/, int
fvHitPointIndexes.push_back(th.iMC);
}
- if (fVerbose >= 1) cout << "ReadEvent: mvd and sts are saved." << endl;
+ if (fVerbose >= 2) cout << "ReadEvent: mvd and sts are saved." << endl;
// ----- Send data from IODataManager to L1Algo --------------------------------------------------------------------
if (1 == fSTAPDataMode) { WriteSTAPAlgoInputData(nCalls); }
diff --git a/reco/L1/L1Algo/L1Constants.h b/reco/L1/L1Algo/L1Constants.h
index 6d8e3cbeed..c11681ad6a 100644
--- a/reco/L1/L1Algo/L1Constants.h
+++ b/reco/L1/L1Algo/L1Constants.h
@@ -72,8 +72,15 @@ namespace L1Constants
/* Particle masses used for track fit */
constexpr float kMuonMass = 0.10565800f; ///< Muon mass [GeV/c2]
constexpr float kElectronMass = 0.000511f; ///< Electron mass [GeV/c2]
+ constexpr double kSpeedOfLight = 29.9792458; ///< Speed of light [cm/ns]
} // namespace phys
+ /// Math constants
+ namespace math
+ {
+ constexpr double kPi = 3.14159265358979323846; ///< Value of PI, used in ROOT TMath
+ }
+
/// Miscellaneous constants
namespace misc
{
diff --git a/reco/L1/L1Algo/L1Undef.h b/reco/L1/L1Algo/L1Undef.h
new file mode 100644
index 0000000000..adc32adf41
--- /dev/null
+++ b/reco/L1/L1Algo/L1Undef.h
@@ -0,0 +1,70 @@
+/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergey Gorbunov, Sergei Zharko [committer] */
+
+/// \file L1Undef.h
+/// \brief Namespace L1Undef provides undefined values for different types
+/// \since 18.11.2022
+/// \author S.Zharko <s.zharko@gsi.de>
+
+#ifndef L1Undef_h
+#define L1Undef_h 1
+
+#include <limits>
+
+#include "CaSimd.h"
+
+// NOTE: SZh 18.11.2022:
+//
+// A temporary solution. Previously, it was found that "using namespace cbm::algo::ca;" was introduced into a L1NaN.h
+// header. This had an influence on every class in L1Algo directory (each of them has at least one fvec or fscal
+// variable). In future classes inside this directory will be placed in the cbm::algo:;ca namespace, so the following
+// lines will not be needed anymore.
+using cbm::algo::ca::fmask;
+using cbm::algo::ca::fscal;
+using cbm::algo::ca::fvec;
+
+namespace undef
+{
+ constexpr int kI32 = -1;
+ constexpr unsigned kU32 = 4294967295;
+ constexpr float kF = std::numeric_limits<float>::quiet_NaN();
+ constexpr double kD = std::numeric_limits<double>::quiet_NaN();
+ constexpr fscal kFvc = std::numeric_limits<fscal>::quiet_NaN();
+
+ /// Checks whether a variable of a particular type defined
+ template<typename T>
+ bool IsUndefined(T val);
+} // namespace undef
+
+template<>
+inline bool undef::IsUndefined<int>(int val)
+{
+ return val == undef::kI32;
+}
+
+template<>
+inline bool undef::IsUndefined<unsigned>(unsigned val)
+{
+ return val == undef::kU32;
+}
+
+template<>
+inline bool undef::IsUndefined<float>(float val)
+{
+ return std::isnan(val);
+}
+
+template<>
+inline bool undef::IsUndefined<double>(double val)
+{
+ return std::isnan(val);
+}
+
+template<>
+inline bool undef::IsUndefined<cbm::algo::ca::fvec>(cbm::algo::ca::fvec val)
+{
+ return isnan(val).isNotEmpty();
+}
+
+#endif // L1Undef_h
diff --git a/reco/L1/L1LinkDef.h b/reco/L1/L1LinkDef.h
index ad266732d7..3681e47ae7 100644
--- a/reco/L1/L1LinkDef.h
+++ b/reco/L1/L1LinkDef.h
@@ -32,6 +32,8 @@
#pragma link C++ class CbmTrackerInputQaTof + ;
#pragma link C++ class CbmCaInputQaMuch + ;
#pragma link C++ class CbmCaInputQaSts + ;
+#pragma link C++ class CbmCaInputQaTrd + ;
+#pragma link C++ class CbmCaInputQaTof + ;
#pragma link C++ class ca::tools::WindowFinder + ;
#endif
diff --git a/reco/L1/catools/CaToolsLinkKey.h b/reco/L1/catools/CaToolsLinkKey.h
new file mode 100644
index 0000000000..aecc629ec4
--- /dev/null
+++ b/reco/L1/catools/CaToolsLinkKey.h
@@ -0,0 +1,52 @@
+/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// \file CaToolsLinkKey.h
+/// \brief Data structure to represent a MC link in CA tracking MC module
+/// \since 25.11.2022
+/// \author S.Zharko <s.zharko@gsi.de>
+
+#ifndef CaToolsLinkKey_h
+#define CaToolsLinkKey_h 1
+
+#include <boost/functional/hash.hpp>
+
+namespace ca::tools
+{
+ struct LinkKey {
+ /// Constructor from index, event and file of a link
+ /// \param index Index of MC track/point in external data structure
+ /// \param event Index of MC event
+ /// \param file Index of MC file
+ LinkKey(int index, int event, int file) : fIndex(index), fEvent(event), fFile(file) {}
+
+ /// Comparison operator
+ friend bool operator==(const LinkKey& lhs, const LinkKey& rhs)
+ {
+ return lhs.fFile == rhs.fFile && lhs.fEvent == rhs.fEvent && lhs.fIndex == rhs.fIndex;
+ }
+
+ int fIndex = -1; ///< Index of MC point/track in external data structures
+ int fEvent = -1; ///< Index of MC event
+ int fFile = -1; ///< Index of MC file
+ };
+} // namespace ca::tools
+
+namespace std
+{
+ /// A hash specialization for ca::tools::LinkKey objects
+ template<>
+ struct hash<ca::tools::LinkKey> {
+ std::size_t operator()(const ca::tools::LinkKey& key) const
+ {
+ std::size_t res = 0;
+ boost::hash_combine(res, key.fFile);
+ boost::hash_combine(res, key.fEvent);
+ boost::hash_combine(res, key.fIndex);
+ return res;
+ }
+ };
+} // namespace std
+
+#endif // CaToolsLinkKey_h
diff --git a/reco/L1/qa/CbmCaInputQaBase.h b/reco/L1/qa/CbmCaInputQaBase.h
new file mode 100644
index 0000000000..fc00be730e
--- /dev/null
+++ b/reco/L1/qa/CbmCaInputQaBase.h
@@ -0,0 +1,80 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// \file CbmCaInputQaBase.h
+/// \date 13.01.2023
+/// \brief Class providing basic functionality for CA input QA-tasks
+/// \author S.Zharko <s.zharko@lx-pool.gsi.de>
+
+#ifndef CbmCaInputQaBase_h
+#define CbmCaInputQaBase_h 1
+
+#include "Logger.h"
+
+#include <array>
+
+/// Class CbmCaInputQaBase contains common functionality for all of the CA input QA-tasks
+///
+class CbmCaInputQaBase {
+public:
+ /// Default constructor
+ CbmCaInputQaBase() = default;
+
+ /// Default destructor
+ ~CbmCaInputQaBase() = default;
+
+
+ /// Sets maximum allowed deviation of pull distribution width from unity
+ /// \param devThrsh Deviation threshold (1 - devThrsh, 1 + devThrsh)
+ void SetPullWidthDeviationThrsh(double devThrsh) { fPullWidthThrsh = devThrsh; }
+
+ /// Sets maximum allowed deviation of residual mean from zero (-devThrsh * rms, +devThrsh * rms)
+ /// \param devThrsh Deviation threshold in sigmas of residual distribution
+ void SetResidualMeanDeviationThrsh(double devThrsh) { fResMeanThrsh = devThrsh; }
+
+ /// Sets hit efficiency analysis parameters: efficiency fit range
+ /// \param thrsh Threshold for efficiency in a given range
+ /// \param loRange Lower bound of the integration range
+ /// \param upRange Upper limit of the integration range
+ void SetEfficiencyThrsh(double thrsh, double loRange, double upRange)
+ {
+ LOG(info) << "Call: SetEfficiencyThrsh";
+ fEffThrsh = thrsh;
+ fEffRange[0] = loRange;
+ fEffRange[1] = upRange;
+ LOG(info) << fEffRange[0] << ' ' << fEffRange[1];
+ }
+
+ /// Sets maximum allowed difference between z-position of hit and station
+ /// \param dz Difference between station and hit position z components [cm]
+ void SetMaxDiffZStationHit(double dz) { fMaxDiffZStHit = dz; }
+
+ /// Sets lower momentum threshold
+ /// \param mom Minimum absolute value of momentum
+ void SetMinMomentum(double mom) { fMinMomentum = mom; }
+
+
+protected:
+ // *********************
+ // ** Data structures **
+ // *********************
+
+ /// Structure to store fit residuals result
+ struct ResidualFitResult {
+ double mean = 0; ///< mean of the distribution
+ double lo = 0; ///< lower limit for the mean
+ double hi = 0; ///< higher limit for the mean
+ };
+
+ // ----- QA constants (thresholds, ranges etc.)
+ double fResMeanThrsh = .5; ///< Maximum allowed deviation of residual mean from zero [sigma]
+ double fPullWidthThrsh = .1; ///< Maximum allowed deviation of pull width from unity
+ double fEffThrsh = .5; ///< Threshold for hit efficiency in the selected range
+ double fMaxDiffZStHit = 1.; ///< Maximum allowed difference between z-position of hit and station [cm]
+ double fMinMomentum = .05; ///< Minimum momentum of particle [GeV/c]
+
+ std::array<double, 2> fEffRange = {0., 100.}; ///< Range for hit efficiency approximation
+};
+
+#endif // CbmCaInputQaBase_h
diff --git a/reco/L1/qa/CbmCaInputQaMuch.cxx b/reco/L1/qa/CbmCaInputQaMuch.cxx
new file mode 100644
index 0000000000..ed8d3d823c
--- /dev/null
+++ b/reco/L1/qa/CbmCaInputQaMuch.cxx
@@ -0,0 +1,1020 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// \file CbmCaInputQaMuch.cxx
+/// \date 13.01.2023
+/// \brief QA-task for CA tracking input from MuCh detector (implementation)
+/// \author S.Zharko <s.zharko@gsi.de>
+
+#include "CbmCaInputQaMuch.h"
+
+#include "CbmAddress.h"
+#include "CbmMCDataArray.h"
+#include "CbmMCEventList.h"
+#include "CbmMCTrack.h"
+#include "CbmMatch.h"
+#include "CbmMuchPixelHit.h"
+#include "CbmMuchPoint.h"
+#include "CbmMuchTrackingInterface.h"
+#include "CbmQaCanvas.h"
+#include "CbmQaEff.h"
+#include "CbmTimeSlice.h"
+
+#include "FairRootManager.h"
+#include "Logger.h"
+
+#include "TBox.h"
+#include "TClonesArray.h"
+#include "TEllipse.h"
+#include "TF1.h"
+#include "TFormula.h"
+#include "TGraphAsymmErrors.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TMath.h"
+#include "TStyle.h"
+
+#include <algorithm>
+#include <fstream>
+#include <iomanip>
+#include <numeric>
+
+#include "L1Constants.h"
+
+ClassImp(CbmCaInputQaMuch);
+
+namespace phys = L1Constants::phys; // from L1Constants.h
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+CbmCaInputQaMuch::CbmCaInputQaMuch(int verbose, bool isMCUsed)
+ : CbmQaTask("CbmCaInputQaMuch", "camuch", verbose, isMCUsed)
+{
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+bool CbmCaInputQaMuch::Check()
+{
+ bool res = true;
+
+ int nSt = fpDetInterface->GetNtrackingStations();
+
+
+ // **************************************************************
+ // ** Basic checks, available both for real and simulated data **
+ // **************************************************************
+
+ // ----- Checks for mismatches in the ordering of the stations
+ //
+ std::vector<double> vStationPos(nSt, 0.);
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ vStationPos[iSt] = fpDetInterface->GetZ(iSt);
+ }
+
+ if (!std::is_sorted(vStationPos.cbegin(), vStationPos.cend(), [](int l, int r) { return l <= r; })) {
+ if (fVerbose > 0) {
+ LOG(error) << fName << ": stations are ordered improperly along the beam axis:";
+ for (auto zPos : vStationPos) {
+ LOG(error) << "\t- " << zPos;
+ }
+ }
+ res = false;
+ }
+
+ // ----- Checks for mismatch between station and hit z positions
+ // The purpose of this block is to be ensured, that hits belong to the correct tracking station. For each tracking
+ // station a unified position z-coordinate is defined, which generally differs from the corresponding positions of
+ // reconstructed hits. This is due to non-regular position along the beam axis for each detector sensor. Thus, the
+ // positions of the hits along z-axis are distributed relatively to the position of the station in some range.
+ // If hits goes out this range, it can signal about a mismatch between hits and geometry. For limits of the range
+ // one should select large enough values to cover the difference described above and in the same time small enough
+ // to avoid overlaps with the neighboring stations.
+ // For each station, a distribution of z_{hit} - z_{st} is integrated over the defined range and scaled by the
+ // total number of entries to the distribution. If this value is smaller then unity, then some of the hits belong to
+ // another station.
+ //
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ int nHits = fvph_hit_station_delta_z[iSt]->GetEntries();
+ if (!nHits) {
+ LOG_IF(error, fVerbose > 0) << fName << ": station " << iSt << " does not have hits";
+ res = false;
+ continue;
+ }
+ int iBinMin = fvph_hit_station_delta_z[iSt]->FindBin(-fMaxDiffZStHit);
+ int iBinMax = fvph_hit_station_delta_z[iSt]->FindBin(+fMaxDiffZStHit);
+
+ if (fvph_hit_station_delta_z[iSt]->Integral(iBinMin, iBinMax) < nHits) {
+ LOG_IF(error, fVerbose > 0) << fName << ": station " << iSt << " has mismatches in hit z-positions";
+ res = false;
+ }
+ }
+
+
+ // *******************************************************
+ // ** Additional checks, if MC information is available **
+ // *******************************************************
+
+ if (IsMCUsed()) {
+ // ----- Check efficiencies
+ // Error is raised, if any station has integrated efficiency lower then a selected threshold.
+ // Integrated efficiency is estimated with fitting the efficiency curve for hit distances (r) with a uniform
+ // distribution in the range from kEffRangeMin to kEffRangeMax, where the efficiency is assigned to be constant
+ //
+ // Fit efficiency curves
+ LOG(info) << "-- Hit efficiency integrated over hit distance from station center";
+ LOG(info) << "\tintegration range: [" << fEffRange[0] << ", " << fEffRange[1] << "] cm";
+
+ auto* pEffTable = MakeTable("eff_table", "Efficiency table", nSt, 2);
+ pEffTable->SetNamesOfCols({"Station ID", "Efficiency"});
+ pEffTable->SetColWidth(20);
+
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ auto pEff = std::unique_ptr<CbmQaEff>(fvpe_reco_eff_vs_r[iSt]->Integral(fEffRange[0], fEffRange[1]));
+ double eff = pEff->GetEfficiency(1);
+ pEffTable->SetCell(iSt, 0, iSt);
+ pEffTable->SetCell(iSt, 1, eff);
+ res = CheckRange("Hit finder efficiency in station " + std::to_string(iSt), eff, fEffThrsh, 1.000);
+ }
+
+ LOG(info) << '\n' << pEffTable->ToString(3);
+
+ // ----- Checks for residuals
+ // Check hits for potential biases from central values
+
+ // Function to fit a residual distribution, returns a structure
+ auto FitResidualsAndGetMean = [&](TH1* pHist) {
+ auto* pFit = (TF1*) gROOT->FindObject("gausn");
+ LOG_IF(fatal, !pFit) << fName << ": residuals fit function is not found";
+ double statMean = pHist->GetMean();
+ double statSigm = pHist->GetStdDev();
+ pFit->SetParameters(100., statMean, statSigm);
+ pHist->Fit(pFit, "MQ");
+ pHist->Fit(pFit, "MQ");
+ pHist->Fit(pFit, "MQ");
+ // NOTE: Several fit procedures are made to avoid empty fit results
+ std::array<double, 3> result;
+ result[0] = pFit->GetParameter(1);
+ result[1] = -pFit->GetParameter(2) * fResMeanThrsh;
+ result[2] = +pFit->GetParameter(2) * fResMeanThrsh;
+ return result;
+ };
+
+ auto* pResidualsTable = MakeTable("residuals_mean", "Residual mean values in different stations", nSt, 4);
+ pResidualsTable->SetNamesOfCols({"Station ID", "Residual(x) [cm]", "Residual(y) [cm]", "Residual(t) [ns]"});
+ pResidualsTable->SetColWidth(20);
+
+ // Fill residuals fit results and check boundaries
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ auto fitX = FitResidualsAndGetMean(fvph_res_x[iSt]);
+ auto fitY = FitResidualsAndGetMean(fvph_res_y[iSt]);
+ auto fitT = FitResidualsAndGetMean(fvph_res_t[iSt]);
+ res = CheckRange("Mean of x residual [cm] in station " + std::to_string(iSt), fitX[0], fitX[1], fitX[2]) && res;
+ res = CheckRange("Mean of y residual [cm] in station " + std::to_string(iSt), fitY[0], fitY[1], fitY[2]) && res;
+ res = CheckRange("Mean of t residual [ns] in station " + std::to_string(iSt), fitT[0], fitT[1], fitT[2]) && res;
+ pResidualsTable->SetCell(iSt, 0, iSt);
+ pResidualsTable->SetCell(iSt, 1, fitX[0]);
+ pResidualsTable->SetCell(iSt, 2, fitX[1]);
+ pResidualsTable->SetCell(iSt, 3, fitX[2]);
+ }
+
+ LOG(info) << '\n' << pResidualsTable->ToString(8);
+
+ // ----- Checks for pulls
+ // For the hit pull is defined as a ration of the difference between hit and MC-point position or time component
+ // values and an error of the hit position (or time) component, calculated in the same z-positions. In ideal case,
+ // when the resolutions are well determined for detecting stations, the pull distributions should have a RMS equal
+ // to unity. Here we allow a RMS of the pull distributions to be varied in a predefined range. If the RMS runs out
+ // this range, QA task fails.
+ // TODO: Add checks for center
+
+ // Fit pull distributions
+
+ // ** Kaniadakis Gaussian distribution, gives smaller chi2 / NDF
+ //if (!gROOT->FindObject("Expk")) {
+ // new TFormula("Expk", "TMath::Power(TMath::Sqrt(1 + x[1] * x[1] * x[0] * x[0]) + x[0] * x[1], 1./x[1])");
+ //}
+ //TF1* pPullFit = new TF1("pullFitGausn", "[0] * Expk(-[2] * (x - [1]) * (x - [1]), [3])", -10., 10.);
+
+ auto FitPullsAndGetSigma = [&](TH1* pHist) {
+ auto* pFit = (TF1*) gROOT->FindObject("gausn");
+ LOG_IF(fatal, !pFit) << fName << ": pulls fit function is not found";
+ pFit->SetParameters(100, 0.001, 1.000);
+ pHist->Fit(pFit, "Q");
+ return pFit->GetParameter(2);
+ };
+
+ auto* pPullsTable = MakeTable("pulls_rms", "Pulls std. dev. values in different stations", nSt, 4);
+ pPullsTable->SetNamesOfCols({"Station ID", "Pull(x) sigm", "Pull(y) sigm", "Pull(z) sigm"});
+ pPullsTable->SetColWidth(20);
+
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ double rmsPullX = FitPullsAndGetSigma(fvph_pull_x[iSt]);
+ double rmsPullY = FitPullsAndGetSigma(fvph_pull_y[iSt]);
+ double rmsPullT = FitPullsAndGetSigma(fvph_pull_t[iSt]);
+
+ double rmsPullHi = 1. + fPullWidthThrsh;
+ double rmsPullLo = 1. - fPullWidthThrsh;
+
+ res = CheckRange("Rms for x pull in station " + std::to_string(iSt), rmsPullX, rmsPullLo, rmsPullHi) && res;
+ res = CheckRange("Rms for y pull in station " + std::to_string(iSt), rmsPullY, rmsPullLo, rmsPullHi) && res;
+ res = CheckRange("Rms for t pull in station " + std::to_string(iSt), rmsPullT, rmsPullLo, rmsPullHi) && res;
+ pPullsTable->SetCell(iSt, 0, iSt);
+ pPullsTable->SetCell(iSt, 1, rmsPullX);
+ pPullsTable->SetCell(iSt, 2, rmsPullY);
+ pPullsTable->SetCell(iSt, 3, rmsPullT);
+ }
+
+ LOG(info) << '\n' << pPullsTable->ToString(3);
+ }
+
+ return res;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void CbmCaInputQaMuch::DeInit()
+{
+ // Vectors with pointers to histograms
+ fvph_hit_ypos_vs_xpos.clear();
+ fvph_hit_xpos_vs_zpos.clear();
+ fvph_hit_ypos_vs_zpos.clear();
+
+ fvph_hit_station_delta_z.clear();
+
+ fvph_hit_dx.clear();
+ fvph_hit_dy.clear();
+ fvph_hit_dt.clear();
+
+ fvph_n_points_per_hit.clear();
+
+ fvph_point_ypos_vs_xpos.clear();
+ fvph_point_xpos_vs_zpos.clear();
+ fvph_point_ypos_vs_zpos.clear();
+
+ fvph_point_hit_delta_z.clear();
+
+ fvph_res_x.clear();
+ fvph_res_y.clear();
+ fvph_res_t.clear();
+
+ fvph_pull_x.clear();
+ fvph_pull_y.clear();
+ fvph_pull_t.clear();
+
+ fvph_res_x_vs_x.clear();
+ fvph_res_y_vs_y.clear();
+ fvph_res_t_vs_t.clear();
+
+ fvph_pull_x_vs_x.clear();
+ fvph_pull_y_vs_y.clear();
+ fvph_pull_t_vs_t.clear();
+
+ fvpe_reco_eff_vs_xy.clear();
+ fvpe_reco_eff_vs_r.clear();
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void CbmCaInputQaMuch::FillHistograms()
+{
+ int nSt = fpDetInterface->GetNtrackingStations();
+ int nHits = fpHits->GetEntriesFast();
+ int nMCevents = (IsMCUsed()) ? fpMCEventList->GetNofEvents() : -1;
+
+ std::vector<std::vector<int>> vNofHitsPerPoint; // Number of hits per MC point in different MC events
+
+ if (IsMCUsed()) {
+ vNofHitsPerPoint.resize(nMCevents);
+ for (int iE = 0; iE < nMCevents; ++iE) {
+ int iFile = fpMCEventList->GetFileIdByIndex(iE);
+ int iEvent = fpMCEventList->GetEventIdByIndex(iE);
+ int nPoints = fpMCPoints->Size(iFile, iEvent);
+ vNofHitsPerPoint[iE].resize(nPoints, 0);
+ }
+ }
+
+ for (int iH = 0; iH < nHits; ++iH) {
+ const auto* pHit = dynamic_cast<const CbmMuchPixelHit*>(fpHits->At(iH));
+ LOG_IF(fatal, !pHit) << fName << ": hit with iH = " << iH << " is not an CbmMuchPixelHit (dynamic cast failed)";
+
+ // *************************
+ // ** Reconstructed hit QA
+
+ // Hit station geometry info
+ int iSt = fpDetInterface->GetTrackingStationIndex(pHit);
+ LOG_IF(fatal, iSt < 0 || iSt >= nSt) << fName << ": index of station (" << iSt << ") is out of range "
+ << "for hit with id = " << iH;
+
+ // Hit position
+ double xHit = pHit->GetX();
+ double yHit = pHit->GetY();
+ double zHit = pHit->GetZ();
+
+ double dxHit = pHit->GetDx();
+ double dyHit = pHit->GetDy();
+ //double dxyHit = pHit->GetDxy();
+
+ // Hit time
+ double tHit = pHit->GetTime();
+ double dtHit = pHit->GetTimeError();
+
+ fvph_hit_ypos_vs_xpos[iSt]->Fill(xHit, yHit);
+ fvph_hit_xpos_vs_zpos[iSt]->Fill(zHit, xHit);
+ fvph_hit_ypos_vs_zpos[iSt]->Fill(zHit, yHit);
+
+ fvph_hit_station_delta_z[iSt]->Fill(zHit - fpDetInterface->GetZ(iSt));
+
+ fvph_hit_dx[iSt]->Fill(dxHit);
+ fvph_hit_dy[iSt]->Fill(dyHit);
+ fvph_hit_dt[iSt]->Fill(dtHit);
+
+ fvph_hit_ypos_vs_xpos[nSt]->Fill(xHit, yHit);
+ fvph_hit_xpos_vs_zpos[nSt]->Fill(zHit, xHit);
+ fvph_hit_ypos_vs_zpos[nSt]->Fill(zHit, yHit);
+ fvph_hit_dx[nSt]->Fill(dxHit);
+ fvph_hit_dy[nSt]->Fill(dyHit);
+ fvph_hit_dt[nSt]->Fill(dtHit);
+
+
+ // **********************
+ // ** MC information QA
+
+ if (IsMCUsed()) {
+ CbmMatch* pHitMatch = dynamic_cast<CbmMatch*>(fpHitMatches->At(iH));
+ LOG_IF(fatal, !pHitMatch) << fName << ": match object not found for hit ID " << iH;
+
+ // Evaluate number of hits per one MC point
+ int nMCpoints = 0; // Number of MC points for this hit
+ for (int iLink = 0; iLink < pHitMatch->GetNofLinks(); ++iLink) {
+ const auto& link = pHitMatch->GetLink(iLink);
+
+ int iP = link.GetIndex(); // Index of MC point
+
+ // Skip noisy links
+ if (iP < 0) { continue; }
+
+ ++nMCpoints;
+
+ int iE = fpMCEventList->GetEventIndex(link);
+
+ LOG_IF(fatal, iE < 0 || iE >= nMCevents) << fName << ": id of MC event is out of range (hit id = " << iH
+ << ", link id = " << iLink << ", event id = " << iE << ')';
+
+ LOG_IF(fatal, iP >= (int) vNofHitsPerPoint[iE].size())
+ << fName << ": index of MC point is out of range (hit id = " << iH << ", link id = " << iLink
+ << ", event id = " << iE << ", point id = " << iP << ')';
+ vNofHitsPerPoint[iE][iP]++;
+ }
+
+ fvph_n_points_per_hit[iSt]->Fill(nMCpoints);
+ fvph_n_points_per_hit[nSt]->Fill(nMCpoints);
+
+ if (nMCpoints != 1) { continue; } // ??
+
+ // The best link to in the match (probably, the cut on nMCpoints is meaningless)
+ const auto& bestPointLink = pHitMatch->GetMatchedLink();
+
+ // Skip noisy links
+ if (bestPointLink.GetIndex() < 0) { continue; }
+
+ // Point matched by the best link
+ const auto* pMCPoint = dynamic_cast<const CbmMuchPoint*>(fpMCPoints->Get(bestPointLink));
+ LOG_IF(fatal, !pMCPoint) << fName << ": MC point object does not exist for hit " << iH;
+
+ // MC track for this point
+ CbmLink bestTrackLink = bestPointLink;
+ bestTrackLink.SetIndex(pMCPoint->GetTrackID());
+ const auto* pMCTrack = dynamic_cast<const CbmMCTrack*>(fpMCTracks->Get(bestTrackLink));
+ LOG_IF(fatal, !pMCTrack) << fName << ": MC track object does not exist for hit " << iH << " and link: ";
+
+ double t0MC = fpMCEventList->GetEventTime(bestPointLink);
+ LOG_IF(fatal, t0MC < 0) << fName << ": MC time zero is lower then 0 ns: " << t0MC;
+
+
+ // ----- MC point properties
+ //
+ double mass = pMCTrack->GetMass();
+ //double charge = pMCTrack->GetCharge();
+ //double pdg = pMCTrack->GetPdgCode();
+
+ // Entrance position and time
+ double xMC = pMCPoint->GetXIn();
+ double yMC = pMCPoint->GetYIn();
+ double zMC = pMCPoint->GetZIn();
+ double tMC = pMCPoint->GetTime() + t0MC;
+
+ // MC point entrance momenta
+ double pxMC = pMCPoint->GetPx();
+ double pyMC = pMCPoint->GetPy();
+ double pzMC = pMCPoint->GetPz();
+ double pMC = sqrt(pxMC * pxMC + pyMC * pyMC + pzMC * pzMC);
+
+ // MC point exit momenta
+ // double pxMCo = pMCPoint->GetPxOut();
+ // double pyMCo = pMCPoint->GetPyOut();
+ // double pzMCo = pMCPoint->GetPzOut();
+ // double pMCo = sqrt(pxMCo * pxMCo + pyMCo * pyMCo + pzMCo * pzMCo);
+
+ // Position and time shifted to z-coordinate of the hit
+ double shiftZ = zHit - zMC; // Difference between hit and point z positions
+ double xMCs = xMC + shiftZ * pxMC / pzMC;
+ double yMCs = yMC + shiftZ * pyMC / pzMC;
+ double tMCs = tMC + shiftZ / (pzMC * phys::kSpeedOfLight) * sqrt(mass * mass + pMC * pMC);
+
+ // Residuals
+ double xRes = xHit - xMCs;
+ double yRes = yHit - yMCs;
+ double tRes = tHit - tMCs;
+
+ // ------ Cuts
+
+ if (std::fabs(pMCPoint->GetPzOut()) < fMinMomentum) { continue; } // CUT ON MINIMUM MOMENTUM
+ //if (pMCo < cuts::kMinP) { continue; } // Momentum threshold
+
+
+ fvph_point_ypos_vs_xpos[iSt]->Fill(xMC, yMC);
+ fvph_point_xpos_vs_zpos[iSt]->Fill(zMC, xMC);
+ fvph_point_ypos_vs_zpos[iSt]->Fill(zMC, yMC);
+
+ fvph_point_hit_delta_z[iSt]->Fill(shiftZ);
+
+ fvph_res_x[iSt]->Fill(xRes);
+ fvph_res_y[iSt]->Fill(yRes);
+ fvph_res_t[iSt]->Fill(tRes);
+
+ fvph_pull_x[iSt]->Fill(xRes / dxHit);
+ fvph_pull_y[iSt]->Fill(yRes / dyHit);
+ fvph_pull_t[iSt]->Fill(tRes / dtHit);
+
+ fvph_res_x_vs_x[iSt]->Fill(xHit, xRes);
+ fvph_res_y_vs_y[iSt]->Fill(yHit, yRes);
+ fvph_res_t_vs_t[iSt]->Fill(tHit, tRes);
+
+ fvph_pull_x_vs_x[iSt]->Fill(xHit, xRes / dxHit);
+ fvph_pull_y_vs_y[iSt]->Fill(yHit, yRes / dyHit);
+ fvph_pull_t_vs_t[iSt]->Fill(tHit, tRes / dtHit);
+
+ fvph_point_ypos_vs_xpos[nSt]->Fill(xMC, yMC);
+ fvph_point_xpos_vs_zpos[nSt]->Fill(zMC, xMC);
+ fvph_point_ypos_vs_zpos[nSt]->Fill(zMC, yMC);
+
+ fvph_point_hit_delta_z[nSt]->Fill(shiftZ);
+
+ fvph_res_x[nSt]->Fill(xRes);
+ fvph_res_y[nSt]->Fill(yRes);
+ fvph_res_t[nSt]->Fill(tRes);
+
+ fvph_pull_x[nSt]->Fill(xRes / dxHit);
+ fvph_pull_y[nSt]->Fill(yRes / dyHit);
+ fvph_pull_t[nSt]->Fill(tRes / dtHit);
+
+ fvph_res_x_vs_x[nSt]->Fill(xHit, xRes);
+ fvph_res_y_vs_y[nSt]->Fill(yHit, yRes);
+ fvph_res_t_vs_t[nSt]->Fill(tHit, tRes);
+
+ fvph_pull_x_vs_x[nSt]->Fill(xHit, xRes / dxHit);
+ fvph_pull_y_vs_y[nSt]->Fill(yHit, yRes / dyHit);
+ fvph_pull_t_vs_t[nSt]->Fill(tHit, tRes / dtHit);
+ }
+ } // loop over hits: end
+
+ // Fill hit reconstruction efficiencies
+ if (IsMCUsed()) {
+ for (int iE = 0; iE < nMCevents; ++iE) {
+ int iFile = fpMCEventList->GetFileIdByIndex(iE);
+ int iEvent = fpMCEventList->GetEventIdByIndex(iE);
+ int nPoints = fpMCPoints->Size(iFile, iEvent);
+
+ for (int iP = 0; iP < nPoints; ++iP) {
+ const auto* pMCPoint = dynamic_cast<const CbmMuchPoint*>(fpMCPoints->Get(iFile, iEvent, iP));
+ LOG_IF(fatal, !pMCPoint) << fName << ": MC point does not exist for iFile = " << iFile
+ << ", iEvent = " << iEvent << ", iP = " << iP;
+ int address = pMCPoint->GetDetectorID();
+ int iSt = fpDetInterface->GetTrackingStationIndex(address);
+ LOG_IF(fatal, iSt < 0 || iSt >= nSt)
+ << fName << ": MC point for FEI = " << iFile << ", " << iEvent << ", " << iP << " and address " << address
+ << " has wrong station index: iSt = " << iSt;
+
+ double xMC = pMCPoint->GetXIn();
+ double yMC = pMCPoint->GetYIn();
+ double rMC = sqrt(xMC * xMC + yMC * yMC);
+
+ // Conditions under which point is accounted as reconstructed: point
+ bool ifPointHasHits = (vNofHitsPerPoint[iE][iP] > 0);
+
+ fvpe_reco_eff_vs_xy[iSt]->Fill(ifPointHasHits, xMC, yMC);
+ fvpe_reco_eff_vs_xy[nSt]->Fill(ifPointHasHits, xMC, yMC);
+
+ fvpe_reco_eff_vs_r[iSt]->Fill(ifPointHasHits, rMC);
+ fvpe_reco_eff_vs_r[nSt]->Fill(ifPointHasHits, rMC);
+
+ } // loop over MC-points: end
+ } // loop over MC-events: end
+ } // MC is used: end
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+InitStatus CbmCaInputQaMuch::InitDataBranches()
+{
+ // STS tracking detector interface
+ fpDetInterface = CbmMuchTrackingInterface::Instance();
+
+ LOG_IF(fatal, !fpDetInterface) << "\033[1;31m" << fName << ": tracking detector interface is undefined\033[0m";
+
+ // FairRootManager
+ auto* pFairRootManager = FairRootManager::Instance();
+ LOG_IF(fatal, !pFairRootManager) << "\033[1;31m" << fName << ": FairRootManager instance is a null pointer\033[0m";
+
+ // Time-slice
+ fpTimeSlice = dynamic_cast<CbmTimeSlice*>(pFairRootManager->GetObject("TimeSlice."));
+ LOG_IF(fatal, !fpTimeSlice) << "\033[1;31m" << fName << ": time-slice branch is not found\033[0m";
+
+ // ----- Reconstructed data-branches initialization
+
+ // Hits container
+ fpHits = dynamic_cast<TClonesArray*>(pFairRootManager->GetObject("MuchPixelHit"));
+ LOG_IF(fatal, !fpHits) << "\033[1;31m" << fName << ": container of reconstructed hits in MuCh is not found\033[0m";
+
+ // ----- MC-information branches initialization
+ if (IsMCUsed()) {
+ // MC manager
+ fpMCDataManager = dynamic_cast<CbmMCDataManager*>(pFairRootManager->GetObject("MCDataManager"));
+ LOG_IF(fatal, !fpMCDataManager) << "\033[1;31m" << fName << ": MC data manager branch is not found\033[0m";
+
+ // MC event list
+ fpMCEventList = dynamic_cast<CbmMCEventList*>(pFairRootManager->GetObject("MCEventList."));
+ LOG_IF(fatal, !fpMCEventList) << "\033[1;31m" << fName << ": MC event list branch is not found\033[0m";
+
+ // MC tracks
+ fpMCTracks = fpMCDataManager->InitBranch("MCTrack");
+ LOG_IF(fatal, !fpMCTracks) << "\033[1;31m" << fName << ": MC track branch is not found\033[0m";
+
+ // MC points
+ fpMCPoints = fpMCDataManager->InitBranch("MuchPoint");
+ LOG_IF(fatal, !fpMCTracks) << "\033[1;31m" << fName << ": MC point branch is not found for MuCh\033[0m";
+
+ // Hit matches
+ fpHitMatches = dynamic_cast<TClonesArray*>(pFairRootManager->GetObject("MuchPixelHitMatch"));
+ LOG_IF(fatal, !fpHitMatches) << "\033[1;31m]" << fName << ": hit match branch is not found for MuCh\033[0m";
+ }
+
+ return kSUCCESS;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+InitStatus CbmCaInputQaMuch::InitCanvases()
+{
+ gStyle->SetOptFit(1);
+ int nSt = fpDetInterface->GetNtrackingStations();
+
+ // ********************
+ // ** Drawing options
+
+ // Contours
+ constexpr auto contColor = kOrange + 7;
+ constexpr auto contWidth = 2; // Line width
+ constexpr auto contStyle = 2; // Line style
+ constexpr auto contFill = 0; // Fill style
+
+ // *********************************
+ // ** Hit occupancies
+
+ // ** Occupancy in XY plane
+ auto* pc_hit_ypos_vs_xpos = MakeCanvas<CbmQaCanvas>("hit_ypos_vs_xpos", "", 1600, 800);
+ pc_hit_ypos_vs_xpos->Divide2D(nSt);
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ pc_hit_ypos_vs_xpos->cd(iSt + 1);
+ fvph_hit_ypos_vs_xpos[iSt]->DrawCopy("colz", "");
+
+
+ // Reference size of the station saved to the detector interface
+ auto* pInnerCircle = new TEllipse(0., 0., fpDetInterface->GetRmin(iSt));
+ pInnerCircle->SetLineWidth(contWidth);
+ pInnerCircle->SetLineStyle(contStyle);
+ pInnerCircle->SetLineColor(contColor);
+ pInnerCircle->SetFillStyle(contFill);
+ pInnerCircle->Draw("SAME");
+
+ auto* pOuterCircle = new TEllipse(0., 0., fpDetInterface->GetRmax(iSt));
+ pOuterCircle->SetLineWidth(contWidth);
+ pOuterCircle->SetLineStyle(contStyle);
+ pOuterCircle->SetLineColor(contColor);
+ pOuterCircle->SetFillStyle(contFill);
+ pOuterCircle->Draw("SAME");
+
+ // Square boarders of station, which are used for the field approximation
+ double stXmin = +fpDetInterface->GetXmax(iSt);
+ double stXmax = -fpDetInterface->GetXmax(iSt);
+ double stYmin = -fpDetInterface->GetYmax(iSt);
+ double stYmax = +fpDetInterface->GetYmax(iSt);
+
+ auto* pBox = new TBox(stXmin, stYmin, stXmax, stYmax);
+ pBox->SetLineWidth(contWidth);
+ pBox->SetLineStyle(contStyle);
+ pBox->SetLineColor(contColor);
+ pBox->SetFillStyle(contFill);
+ pBox->Draw("SAME");
+ }
+
+ // ----- Occupancy projections
+ auto* pc_hit_xpos = MakeCanvas<CbmQaCanvas>("hit_xpos", "", 1400, 800);
+ pc_hit_xpos->Divide2D(nSt);
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ pc_hit_xpos->cd(iSt + 1);
+ auto* phProj = fvph_hit_ypos_vs_xpos[iSt]->ProjectionY();
+ phProj->DrawCopy();
+ }
+
+ auto* pc_hit_ypos = MakeCanvas<CbmQaCanvas>("hit_ypos", "", 1400, 800);
+ pc_hit_ypos->Divide2D(nSt);
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ pc_hit_ypos->cd(iSt + 1);
+ auto* phProj = fvph_hit_ypos_vs_xpos[iSt]->ProjectionX();
+ phProj->DrawCopy();
+ }
+
+
+ // ** Occupancy in XZ plane
+ MakeCanvas<CbmQaCanvas>("hit_xpos_vs_zpos", "", 600, 400);
+ fvph_hit_xpos_vs_zpos[nSt]->DrawCopy("colz", "");
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ // Station positions in detector IFS
+ double stZmin = fpDetInterface->GetZ(iSt) - 0.5 * fpDetInterface->GetThickness(iSt);
+ double stZmax = fpDetInterface->GetZ(iSt) + 0.5 * fpDetInterface->GetThickness(iSt);
+ double stHmin = -fpDetInterface->GetRmax(iSt);
+ double stHmax = +fpDetInterface->GetRmax(iSt);
+
+ auto* pBox = new TBox(stZmin, stHmin, stZmax, stHmax);
+ pBox->SetLineWidth(contWidth);
+ pBox->SetLineStyle(contStyle);
+ pBox->SetLineColor(contColor);
+ pBox->SetFillStyle(contFill);
+ pBox->Draw("SAME");
+ }
+
+ // ** Occupancy in YZ plane
+ MakeCanvas<CbmQaCanvas>("hit_ypos_vs_zpos", "", 600, 400);
+ fvph_hit_ypos_vs_zpos[nSt]->DrawCopy("colz", "");
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ // Station positions in detector IFS
+ double stZmin = fpDetInterface->GetZ(iSt) - 0.5 * fpDetInterface->GetThickness(iSt);
+ double stZmax = fpDetInterface->GetZ(iSt) + 0.5 * fpDetInterface->GetThickness(iSt);
+ double stHmin = -fpDetInterface->GetRmax(iSt);
+ double stHmax = +fpDetInterface->GetRmax(iSt);
+
+ auto* pBox = new TBox(stZmin, stHmin, stZmax, stHmax);
+ pBox->SetLineWidth(contWidth);
+ pBox->SetLineStyle(contStyle);
+ pBox->SetLineColor(contColor);
+ pBox->SetFillStyle(contFill);
+ pBox->Draw("SAME");
+ }
+
+ // ************
+ // ************
+
+ if (IsMCUsed()) {
+
+ // **********************
+ // ** Point occupancies
+
+ // ** Occupancy in XY plane
+ auto* pc_point_ypos_vs_xpos = MakeCanvas<CbmQaCanvas>("point_ypos_vs_xpos", "", 1600, 800);
+ pc_point_ypos_vs_xpos->Divide2D(nSt);
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ pc_point_ypos_vs_xpos->cd(iSt + 1);
+ fvph_point_ypos_vs_xpos[iSt]->DrawCopy("colz", "");
+
+ // Reference size of the station saved to the detector interface
+ auto* pInnerCircle = new TEllipse(0., 0., fpDetInterface->GetRmin(iSt));
+ pInnerCircle->SetLineWidth(contWidth);
+ pInnerCircle->SetLineStyle(contStyle);
+ pInnerCircle->SetLineColor(contColor);
+ pInnerCircle->SetFillStyle(contFill);
+ pInnerCircle->Draw("SAME");
+
+ auto* pOuterCircle = new TEllipse(0., 0., fpDetInterface->GetRmax(iSt));
+ pOuterCircle->SetLineWidth(contWidth);
+ pOuterCircle->SetLineStyle(contStyle);
+ pOuterCircle->SetLineColor(contColor);
+ pOuterCircle->SetFillStyle(contFill);
+ pOuterCircle->Draw("SAME");
+
+ // Square boarders of station, which are used for the field approximation
+ double stXmin = +fpDetInterface->GetXmax(iSt);
+ double stXmax = -fpDetInterface->GetXmax(iSt);
+ double stYmin = -fpDetInterface->GetYmax(iSt);
+ double stYmax = +fpDetInterface->GetYmax(iSt);
+
+ auto* pBox = new TBox(stXmin, stYmin, stXmax, stYmax);
+ pBox->SetLineWidth(contWidth);
+ pBox->SetLineStyle(contStyle);
+ pBox->SetLineColor(contColor);
+ pBox->SetFillStyle(contFill);
+ pBox->Draw("SAME");
+ }
+
+ // ** Occupancy in XZ plane
+ MakeCanvas<CbmQaCanvas>("point_xpos_vs_zpos", "", 600, 400);
+ fvph_point_xpos_vs_zpos[nSt]->SetStats(false);
+ fvph_point_xpos_vs_zpos[nSt]->DrawCopy("colz", "");
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ // Station positions in detector IFS
+ double stZmin = fpDetInterface->GetZ(iSt) - 0.5 * fpDetInterface->GetThickness(iSt);
+ double stZmax = fpDetInterface->GetZ(iSt) + 0.5 * fpDetInterface->GetThickness(iSt);
+ double stHmin = -fpDetInterface->GetRmax(iSt);
+ double stHmax = +fpDetInterface->GetRmax(iSt);
+
+ auto* pBox = new TBox(stZmin, stHmin, stZmax, stHmax);
+ pBox->SetLineWidth(contWidth);
+ pBox->SetLineStyle(contStyle);
+ pBox->SetLineColor(contColor);
+ pBox->SetFillStyle(contFill);
+ pBox->Draw("SAME");
+ }
+
+ // ** Occupancy in YZ plane
+ MakeCanvas<CbmQaCanvas>("point_ypos_vs_zpos", "", 600, 400);
+ fvph_point_ypos_vs_zpos[nSt]->SetStats(false);
+ fvph_point_ypos_vs_zpos[nSt]->DrawCopy("colz", "");
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ // Station positions in detector IFS
+ double stZmin = fpDetInterface->GetZ(iSt) - 0.5 * fpDetInterface->GetThickness(iSt);
+ double stZmax = fpDetInterface->GetZ(iSt) + 0.5 * fpDetInterface->GetThickness(iSt);
+ double stHmin = -fpDetInterface->GetRmax(iSt);
+ double stHmax = +fpDetInterface->GetRmax(iSt);
+
+ auto* pBox = new TBox(stZmin, stHmin, stZmax, stHmax);
+ pBox->SetLineWidth(contWidth);
+ pBox->SetLineStyle(contStyle);
+ pBox->SetLineColor(contColor);
+ pBox->SetFillStyle(contFill);
+ pBox->Draw("SAME");
+ }
+
+ // **************
+ // ** Residuals
+
+ // x-coordinate
+ auto* pc_res_x = MakeCanvas<CbmQaCanvas>("res_x", "Residuals for x coordinate", 1600, 800);
+ pc_res_x->Divide2D(nSt);
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ pc_res_x->cd(iSt + 1);
+ fvph_res_x[iSt]->DrawCopy("", "");
+ }
+
+ // y-coordinate
+ auto* pc_res_y = MakeCanvas<CbmQaCanvas>("res_y", "Residuals for y coordinate", 1600, 800);
+ pc_res_y->Divide2D(nSt);
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ pc_res_y->cd(iSt + 1);
+ fvph_res_y[iSt]->DrawCopy("", "");
+ }
+
+ // time
+ auto* pc_res_t = MakeCanvas<CbmQaCanvas>("res_t", "Residuals for time", 1600, 800);
+ pc_res_t->Divide2D(nSt);
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ pc_res_t->cd(iSt + 1);
+ fvph_res_t[iSt]->DrawCopy("", "");
+ }
+
+
+ // **********
+ // ** Pulls
+
+ // x-coordinate
+ auto* pc_pull_x = MakeCanvas<CbmQaCanvas>("pull_x", "Pulls for x coordinate", 1600, 800);
+ pc_pull_x->Divide2D(nSt);
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ pc_pull_x->cd(iSt + 1);
+ fvph_pull_x[iSt]->DrawCopy("", "");
+ }
+
+ // y-coordinate
+ auto* pc_pull_y = MakeCanvas<CbmQaCanvas>("pull_y", "Pulls for y coordinate", 1600, 800);
+ pc_pull_y->Divide2D(nSt);
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ pc_pull_y->cd(iSt + 1);
+ fvph_pull_y[iSt]->DrawCopy("", "");
+ }
+
+ // time
+ auto* pc_pull_t = MakeCanvas<CbmQaCanvas>("pull_t", "Pulls for time", 1600, 800);
+ pc_pull_t->Divide2D(nSt);
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ pc_pull_t->cd(iSt + 1);
+ fvph_pull_t[iSt]->DrawCopy("", "");
+ }
+
+
+ // ************************************
+ // ** Hit reconstruction efficiencies
+
+ auto* pc_reco_eff_vs_r =
+ MakeCanvas<CbmQaCanvas>("reco_eff_vs_r", "Hit efficiencies wrt distance from center", 1600, 800);
+ pc_reco_eff_vs_r->Divide2D(nSt);
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ pc_reco_eff_vs_r->cd(iSt + 1);
+ fvpe_reco_eff_vs_r[iSt]->Paint("AP");
+ auto* pGr = dynamic_cast<TGraphAsymmErrors*>(fvpe_reco_eff_vs_r[iSt]->GetPaintedGraph());
+ if (!pGr) {
+ LOG(error) << fName << ": unable to get painted graph from efficiency " << fvpe_reco_eff_vs_xy[iSt]->GetName();
+ continue;
+ }
+ pGr->DrawClone("AP");
+ }
+
+ auto* pc_reco_eff_vs_xy = MakeCanvas<CbmQaCanvas>("reco_eff_vs_xy", "Hit efficiencies wrt x and y", 1600, 800);
+ pc_reco_eff_vs_xy->Divide2D(nSt);
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ pc_reco_eff_vs_xy->cd(iSt + 1);
+ fvpe_reco_eff_vs_xy[iSt]->Paint("colz");
+ auto* pH2 = dynamic_cast<TH2F*>(fvpe_reco_eff_vs_xy[iSt]->GetPaintedHistogram());
+ if (!pH2) {
+ LOG(error) << fName << ": unable to get painted histogram from efficiency "
+ << fvpe_reco_eff_vs_xy[iSt]->GetName();
+ continue;
+ }
+ pH2->DrawCopy("colz", "");
+ }
+ }
+
+ return kSUCCESS;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+InitStatus CbmCaInputQaMuch::InitHistograms()
+{
+ int nSt = fpDetInterface->GetNtrackingStations();
+
+ // ----- Histograms initialization
+ fvph_hit_ypos_vs_xpos.resize(nSt + 1, nullptr);
+ fvph_hit_ypos_vs_zpos.resize(nSt + 1, nullptr);
+ fvph_hit_xpos_vs_zpos.resize(nSt + 1, nullptr);
+
+ fvph_hit_station_delta_z.resize(nSt);
+
+ fvph_hit_dx.resize(nSt + 1, nullptr);
+ fvph_hit_dy.resize(nSt + 1, nullptr);
+ fvph_hit_dt.resize(nSt + 1, nullptr);
+
+ for (int iSt = 0; iSt <= nSt; ++iSt) {
+ TString nsuff = (iSt == nSt) ? "" : Form("_st%d", iSt); // Histogram name suffix
+ TString tsuff = (iSt == nSt) ? "" : Form(" in MuCh station %d", iSt); // Histogram title suffix
+ TString sN = "";
+ TString sT = "";
+
+ // Hit occupancy
+ sN = (TString) "hit_ypos_vs_xpos" + nsuff;
+ sT = (TString) "Hit occupancy in xy-Plane" + tsuff + ";x_{hit} [cm];y_{hit} [cm]";
+ fvph_hit_ypos_vs_xpos[iSt] = MakeHisto<TH2F>(sN, sT, kNbins, kRHitX[0], kRHitX[1], kNbins, kRHitY[0], kRHitY[1]);
+
+ sN = (TString) "hit_xpos_vs_zpos" + nsuff;
+ sT = (TString) "Hit occupancy in xz-Plane" + tsuff + ";z_{hit} [cm];x_{hit} [cm]";
+ fvph_hit_xpos_vs_zpos[iSt] = MakeHisto<TH2F>(sN, sT, kNbinsZ, kRHitZ[0], kRHitZ[1], kNbins, kRHitX[0], kRHitX[1]);
+
+ sN = (TString) "hit_ypos_vs_zpos" + nsuff;
+ sT = (TString) "Hit occupancy in yz-plane" + tsuff + ";z_{hit} [cm];y_{hit} [cm]";
+ fvph_hit_ypos_vs_zpos[iSt] = MakeHisto<TH2F>(sN, sT, kNbinsZ, kRHitZ[0], kRHitZ[1], kNbins, kRHitY[0], kRHitY[1]);
+
+ // Hit errors
+ sN = (TString) "hit_dx" + nsuff;
+ sT = (TString) "Hit position error along x-axis" + tsuff + ";dx_{hit} [cm]";
+ fvph_hit_dx[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, kRHitDx[0], kRHitDx[1]);
+
+ sN = (TString) "hit_dy" + nsuff;
+ sT = (TString) "Hit position error along y-axis" + tsuff + ";dy_{hit} [cm]";
+ fvph_hit_dy[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, kRHitDy[0], kRHitDy[1]);
+
+ sN = (TString) "hit_dt" + nsuff;
+ sT = (TString) "Hit time error" + tsuff + ";dt_{hit} [ns]";
+ fvph_hit_dt[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, kRHitDt[0], kRHitDt[1]);
+
+ if (iSt == nSt) { continue; } // Following histograms are defined only for separate station
+
+ sN = (TString) "hit_station_delta_z" + nsuff;
+ sT = (TString) "Different between hit and station z-positions" + tsuff + ";z_{hit} - z_{st} [cm]";
+ fvph_hit_station_delta_z[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, -5., 5);
+
+ } // loop over station index: end
+
+ // ----- Initialize histograms, which are use MC-information
+ if (IsMCUsed()) {
+ // Resize histogram vectors
+ fvph_n_points_per_hit.resize(nSt + 1, nullptr);
+ fvph_point_ypos_vs_xpos.resize(nSt + 1, nullptr);
+ fvph_point_xpos_vs_zpos.resize(nSt + 1, nullptr);
+ fvph_point_ypos_vs_zpos.resize(nSt + 1, nullptr);
+ fvph_point_hit_delta_z.resize(nSt + 1, nullptr);
+ fvph_res_x.resize(nSt + 1, nullptr);
+ fvph_res_y.resize(nSt + 1, nullptr);
+ fvph_res_t.resize(nSt + 1, nullptr);
+ fvph_pull_x.resize(nSt + 1, nullptr);
+ fvph_pull_y.resize(nSt + 1, nullptr);
+ fvph_pull_t.resize(nSt + 1, nullptr);
+ fvph_res_x_vs_x.resize(nSt + 1, nullptr);
+ fvph_res_y_vs_y.resize(nSt + 1, nullptr);
+ fvph_res_t_vs_t.resize(nSt + 1, nullptr);
+ fvph_pull_x_vs_x.resize(nSt + 1, nullptr);
+ fvph_pull_y_vs_y.resize(nSt + 1, nullptr);
+ fvph_pull_t_vs_t.resize(nSt + 1, nullptr);
+ fvpe_reco_eff_vs_xy.resize(nSt + 1, nullptr);
+ fvpe_reco_eff_vs_r.resize(nSt + 1, nullptr);
+
+ for (int iSt = 0; iSt <= nSt; ++iSt) {
+ TString nsuff = (iSt == nSt) ? "" : Form("_st%d", iSt); // Histogram name suffix
+ TString tsuff = (iSt == nSt) ? "" : Form(" in MuCh station %d", iSt); // Histogram title suffix
+ TString sN = "";
+ TString sT = "";
+
+ sN = (TString) "n_points_per_hit" + nsuff;
+ sT = (TString) "Number of points per hit" + tsuff + ";N_{point}/hit";
+ fvph_n_points_per_hit[iSt] = MakeHisto<TH1F>(sN, sT, 10, -0.5, 9.5);
+
+ // Point occupancy
+ sN = (TString) "point_ypos_vs_xpos" + nsuff;
+ sT = (TString) "Point occupancy in XY plane" + tsuff + ";x_{MC} [cm];y_{MC} [cm]";
+ fvph_point_ypos_vs_xpos[iSt] =
+ MakeHisto<TH2F>(sN, sT, kNbins, kRHitX[0], kRHitX[1], kNbins, kRHitY[0], kRHitY[1]);
+
+ sN = (TString) "point_xpos_vs_zpos" + nsuff;
+ sT = (TString) "Point Occupancy in XZ plane" + tsuff + ";z_{MC} [cm];x_{MC} [cm]";
+ fvph_point_xpos_vs_zpos[iSt] =
+ MakeHisto<TH2F>(sN, sT, kNbinsZ, kRHitZ[0], kRHitZ[1], kNbins, kRHitX[0], kRHitX[1]);
+
+ sN = (TString) "point_ypos_vs_zpos" + nsuff;
+ sT = (TString) "Point Occupancy in YZ Plane" + tsuff + ";z_{MC} [cm];y_{MC} [cm]";
+ fvph_point_ypos_vs_zpos[iSt] =
+ MakeHisto<TH2F>(sN, sT, kNbinsZ, kRHitZ[0], kRHitZ[1], kNbins, kRHitY[0], kRHitY[1]);
+
+ // Difference between MC input z and hit z coordinates
+ sN = (TString) "point_hit_delta_z" + nsuff;
+ sT = (TString) "Distance between point and hit along z axis" + tsuff + ";z_{reco} - z_{MC} [cm]";
+ fvph_point_hit_delta_z[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, -0.04, 0.04);
+
+ sN = (TString) "res_x" + nsuff;
+ sT = (TString) "Residuals for x-coordinate" + tsuff + ";x_{reco} - x_{MC} [cm]";
+ fvph_res_x[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, kRResX[0], kRResX[1]);
+
+ sN = (TString) "res_y" + nsuff;
+ sT = (TString) "Residuals for y-coordinate" + tsuff + ";y_{reco} - y_{MC} [cm]";
+ fvph_res_y[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, kRResY[0], kRResY[1]);
+
+ sN = (TString) "res_t" + nsuff;
+ sT = (TString) "Residuals for time" + tsuff + ";t_{reco} - t_{MC} [ns]";
+ fvph_res_t[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, kRResT[0], kRResT[1]);
+
+ sN = (TString) "pull_x" + nsuff;
+ sT = (TString) "Pulls for x-coordinate" + tsuff + ";(x_{reco} - x_{MC}) / #sigma_{x}^{reco}";
+ fvph_pull_x[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, kRPullX[0], kRPullX[1]);
+
+ sN = (TString) "pull_y" + nsuff;
+ sT = (TString) "Pulls for y-coordinate" + tsuff + ";(y_{reco} - y_{MC}) / #sigma_{y}^{reco}";
+ fvph_pull_y[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, kRPullY[0], kRPullY[1]);
+
+ sN = (TString) "pull_t" + nsuff;
+ sT = (TString) "Pulls for time" + tsuff + ";(t_{reco} - t_{MC}) / #sigma_{t}^{reco}";
+ fvph_pull_t[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, kRPullT[0], kRPullT[1]);
+
+ sN = (TString) "res_x_vs_x" + nsuff;
+ sT = (TString) "Residuals for x-coordinate" + tsuff + ";x_{reco} [cm];x_{reco} - x_{MC} [cm]";
+ fvph_res_x_vs_x[iSt] = MakeHisto<TH2F>(sN, sT, kNbins, 0, 0, kNbins, kRResX[0], kRResX[1]);
+
+ sN = (TString) "res_y_vs_y" + nsuff;
+ sT = (TString) "Residuals for y-coordinate" + tsuff + ";y_{reco} [cm];y_{reco} - y_{MC} [cm]";
+ fvph_res_y_vs_y[iSt] = MakeHisto<TH2F>(sN, sT, kNbins, 0, 0, kNbins, kRResY[0], kRResY[1]);
+
+ sN = (TString) "res_t_vs_t" + nsuff;
+ sT = (TString) "Residuals for time" + tsuff + "t_{reco} [ns];t_{reco} - t_{MC} [ns]";
+ fvph_res_t_vs_t[iSt] = MakeHisto<TH2F>(sN, sT, kNbins, 0, 0, kNbins, kRResT[0], kRResT[1]);
+
+ sN = (TString) "pull_x_vs_x" + nsuff;
+ sT = (TString) "Pulls for x-coordinate" + tsuff + "x_{reco} [cm];(x_{reco} - x_{MC}) / #sigma_{x}^{reco}";
+ fvph_pull_x_vs_x[iSt] = MakeHisto<TH2F>(sN, sT, kNbins, 0, 0, kNbins, kRPullX[0], kRPullX[1]);
+
+ sN = (TString) "pull_y_vs_y" + nsuff;
+ sT = (TString) "Pulls for y-coordinate" + tsuff + ";y_{reco} [cm];(y_{reco} - y_{MC}) / #sigma_{y}^{reco}";
+ fvph_pull_y_vs_y[iSt] = MakeHisto<TH2F>(sN, sT, kNbins, 0, 0, kNbins, kRPullY[0], kRPullY[1]);
+
+ sN = (TString) "pull_t_vs_t" + nsuff;
+ sT = (TString) "Pulls for time" + tsuff + ";t_{reco} [ns];(t_{reco} - t_{MC}) / #sigma_{t}^{reco}";
+ fvph_pull_t_vs_t[iSt] = MakeHisto<TH2F>(sN, sT, kNbins, 0, 0, kNbins, kRPullT[0], kRPullT[1]);
+
+
+ sN = (TString) "reco_eff_vs_xy" + nsuff;
+ sT = (TString) "Hit rec. efficiency" + tsuff + ";x_{MC} [cm];y_{MC} [cm];#epsilon";
+ fvpe_reco_eff_vs_xy[iSt] = MakeEfficiency<CbmQaEff>(sN, sT, 100, -50, 50, 100, -50, 50);
+
+ sN = (TString) "reco_eff_vs_r" + nsuff;
+ sT = (TString) "Hit rec. efficiency" + tsuff + ";r_{MC} [cm];#epsilon";
+ fvpe_reco_eff_vs_r[iSt] = MakeEfficiency<CbmQaEff>(sN, sT, 100, 0, 100);
+ }
+ }
+ return kSUCCESS;
+}
diff --git a/reco/L1/qa/CbmCaInputQaMuch.h b/reco/L1/qa/CbmCaInputQaMuch.h
new file mode 100644
index 0000000000..48e65d3853
--- /dev/null
+++ b/reco/L1/qa/CbmCaInputQaMuch.h
@@ -0,0 +1,166 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// \file CbmCaInputQaMuch.h
+/// \date 13.01.2023
+/// \brief QA-task for CA tracking input from MuCh detector (header)
+/// \author S.Zharko <s.zharko@gsi.de>
+
+
+#ifndef CbmCaInputQaMuch_h
+#define CbmCaInputQaMuch_h 1
+
+#include "CbmCaInputQaBase.h"
+#include "CbmMCDataManager.h"
+#include "CbmQaTask.h"
+
+#include "TMath.h"
+
+#include <set>
+#include <unordered_map>
+#include <vector>
+
+class CbmMCEventList;
+class CbmMCDataArray;
+class CbmMCDataManager;
+class CbmTimeSlice;
+class CbmMatch;
+class CbmMuchPixelHit;
+class CbmMuchTrackingInterface;
+class TClonesArray;
+class TH1F;
+class CbmQaEff;
+
+/// A QA-task class, which provides assurance of MuCh hits and geometry
+class CbmCaInputQaMuch : public CbmQaTask, public CbmCaInputQaBase {
+public:
+ /// Constructor from parameters
+ /// \param verbose Verbose level
+ /// \param isMCUsed Flag, whether MC information is available for this task
+ CbmCaInputQaMuch(int verbose, bool isMCUsed);
+
+ ClassDef(CbmCaInputQaMuch, 0);
+
+protected:
+ // ********************************************
+ // ** Virtual method override from CbmQaTask **
+ // ********************************************
+
+ /// Checks results of the QA
+ /// \return Success flag
+ bool Check();
+
+ /// Initializes data branches
+ InitStatus InitDataBranches();
+
+ /// Initializes canvases
+ InitStatus InitCanvases();
+
+ /// Initializes histograms
+ InitStatus InitHistograms();
+
+ /// Fills histograms per event or time-slice
+ void FillHistograms();
+
+ /// De-initializes histograms
+ void DeInit();
+
+private:
+ // *********************
+ // ** Private methods **
+ // *********************
+
+ // *********************
+ // ** Class variables **
+ // *********************
+
+
+ // ** Data branches **
+ CbmMuchTrackingInterface* fpDetInterface = nullptr; ///< Instance of detector interface
+
+ CbmTimeSlice* fpTimeSlice = nullptr; ///< Pointer to current time-slice
+
+ TClonesArray* fpHits = nullptr; ///< Array of hits
+
+ CbmMCDataManager* fpMCDataManager = nullptr; ///< MC data manager
+ CbmMCEventList* fpMCEventList = nullptr; ///< MC event list
+ CbmMCDataArray* fpMCTracks = nullptr; ///< Array of MC tracks
+ CbmMCDataArray* fpMCPoints = nullptr; ///< Array of MC points
+
+ TClonesArray* fpHitMatches = nullptr; ///< Array of hit matches
+
+ // ** Histograms **
+
+ static constexpr double kEffRangeMin = 10.; ///< Lower limit of hit distance for efficiency integration [cm]
+ static constexpr double kEffRangeMax = 30.; ///< Upper limit of hit distance for efficiency integration [cm]
+
+ static constexpr double kRHitX[2] = {-100., 100}; ///< Range for hit x coordinate [cm]
+ static constexpr double kRHitY[2] = {-100., 100}; ///< Range for hit y coordinate [cm]
+ static constexpr double kRHitZ[2] = {145., 335}; ///< Range for hit z coordinate [cm]
+
+ static constexpr int kNbins = 200; ///< General number of bins
+ static constexpr int kNbinsZ = 800; ///< Number of bins for z coordinate axis
+
+ static constexpr double kRHitDx[2] = {-.05, .005}; ///< Range for hit x coordinate [cm]
+ static constexpr double kRHitDy[2] = {-.05, .005}; ///< Range for hit y coordinate [cm]
+ static constexpr double kRHitDt[2] = {-10., 10.}; ///< Range for hit time [ns]
+
+ static constexpr double kRResX[2] = {-.02, .02}; ///< Range for residual of x coordinate [cm]
+ static constexpr double kRResY[2] = {-.1, .1}; ///< Range for residual of y coordinate [cm]
+ static constexpr double kRResT[2] = {-15., 15.}; ///< Range for residual of time [ns]
+
+ static constexpr double kRPullX[2] = {-10., 10.}; ///< Range for pull of x coordinate
+ static constexpr double kRPullY[2] = {-10., 10.}; ///< Range for pull of y coordinate
+ static constexpr double kRPullT[2] = {-5., 5.}; ///< Range for pull of time
+
+
+ // NOTE: the last element of each vector stands for integral distribution over all stations
+
+ // hit occupancy
+ std::vector<TH2F*> fvph_hit_ypos_vs_xpos; ///< hit ypos vs xpos in different stations
+ std::vector<TH2F*> fvph_hit_xpos_vs_zpos; ///< hit xpos vs zpos in different stations
+ std::vector<TH2F*> fvph_hit_ypos_vs_zpos; ///< hit ypos vs zpos in different stations
+
+ // difference between hit and station z positions
+ std::vector<TH1F*> fvph_hit_station_delta_z; ///< Difference between station and hit z positions [cm]
+
+ // hit errors
+ std::vector<TH1F*> fvph_hit_dx;
+ std::vector<TH1F*> fvph_hit_dy;
+ std::vector<TH1F*> fvph_hit_dt;
+
+ // MC points occupancy
+ std::vector<TH1F*> fvph_n_points_per_hit; ///< number of points per one hit in different stations
+
+
+ std::vector<TH2F*> fvph_point_ypos_vs_xpos; ///< point ypos vs xpos in different stations
+ std::vector<TH2F*> fvph_point_xpos_vs_zpos; ///< point xpos vs zpos in different stations
+ std::vector<TH2F*> fvph_point_ypos_vs_zpos; ///< point ypos vs zpos in different stations
+
+ std::vector<TH1F*> fvph_point_hit_delta_z; ///< difference between z of hit and MC point in different stations
+
+ // Residuals
+ std::vector<TH1F*> fvph_res_x; ///< residual for x coordinate in different stations
+ std::vector<TH1F*> fvph_res_y; ///< residual for y coordinate in different stations
+ std::vector<TH1F*> fvph_res_t; ///< residual for t coordinate in different stations
+
+ std::vector<TH2F*> fvph_res_x_vs_x; ///< residual for x coordinate in different stations
+ std::vector<TH2F*> fvph_res_y_vs_y; ///< residual for y coordinate in different stations
+ std::vector<TH2F*> fvph_res_t_vs_t; ///< residual for t coordinate in different stations
+
+ // Pulls
+ std::vector<TH1F*> fvph_pull_x; ///< pull for x coordinate in different stations
+ std::vector<TH1F*> fvph_pull_y; ///< pull for y coordinate in different stations
+ std::vector<TH1F*> fvph_pull_t; ///< pull for t coordinate in different stations
+
+ std::vector<TH2F*> fvph_pull_x_vs_x; ///< pull for x coordinate in different stations
+ std::vector<TH2F*> fvph_pull_y_vs_y; ///< pull for y coordinate in different stations
+ std::vector<TH2F*> fvph_pull_t_vs_t; ///< pull for t coordinate in different stations
+
+ // Hit efficiencies
+ std::vector<CbmQaEff*> fvpe_reco_eff_vs_xy; ///< Efficiency of hit reconstruction vs. x and y coordinates of MC
+ std::vector<CbmQaEff*> fvpe_reco_eff_vs_r; ///< Efficiency of hit reconstruction vs. distance from center
+};
+
+#endif // CbmCaInputQaMuch_h
diff --git a/reco/L1/qa/CbmCaInputQaSts.cxx b/reco/L1/qa/CbmCaInputQaSts.cxx
index 7193fd7f17..a616e002e4 100644
--- a/reco/L1/qa/CbmCaInputQaSts.cxx
+++ b/reco/L1/qa/CbmCaInputQaSts.cxx
@@ -10,13 +10,13 @@
#include "CbmCaInputQaSts.h"
#include "CbmAddress.h"
-#include "CbmCaQaCuts.h"
#include "CbmMCDataArray.h"
#include "CbmMCEventList.h"
#include "CbmMCTrack.h"
#include "CbmMatch.h"
#include "CbmQaCanvas.h"
#include "CbmQaEff.h"
+#include "CbmQaTable.h"
#include "CbmStsCluster.h"
#include "CbmStsHit.h"
#include "CbmStsPoint.h"
@@ -41,63 +41,24 @@
#include <fstream>
#include <iomanip>
#include <numeric>
+#include <tuple>
#include "L1Constants.h"
ClassImp(CbmCaInputQaSts);
-namespace cuts = cbm::ca::qa::cuts; // from CbmCaQaCuts.h
namespace phys = L1Constants::phys; // from L1Constants.h
// ---------------------------------------------------------------------------------------------------------------------
//
-CbmCaInputQaSts::CbmCaInputQaSts(int verbose, bool isMCUsed) : CbmQaTask("CbmCaInputQaSts", verbose, isMCUsed) {}
-
-// ---------------------------------------------------------------------------------------------------------------------
-//
-void CbmCaInputQaSts::AnalyzeHistograms()
+CbmCaInputQaSts::CbmCaInputQaSts(int verbose, bool isMCUsed) : CbmQaTask("CbmCaInputQaSts", "casts", verbose, isMCUsed)
{
- int nSt = fpDetInterface->GetNtrackingStations();
-
- // ----- Fit efficiency curves
- TF1* pFitFn = new TF1("fitfn", "pol0", kEffRangeMin, kEffRangeMax);
-
- for (int iSt = 0; iSt < nSt; ++iSt) {
- pFitFn->SetParameter(0, 0.5);
- fvpe_reco_eff_vs_r[iSt]->Fit(pFitFn);
- }
-
- // ----- Fit pull distributions
- if (!gROOT->FindObject("Expk")) {
- new TFormula("Expk", "TMath::Power(TMath::Sqrt(1 + x[1] * x[1] * x[0] * x[0]) + x[0] * x[1], 1./x[1])");
- }
-
- TF1* pPullFit = new TF1("pullFitGausn", "gausn", -10., 10.);
- //TF1* pPullFit = new TF1("pullFitGausn", "[0] * Expk(-[2] * (x - [1]) * (x - [1]), [3])", -10., 10.);
-
- auto FitPulls = [&](TH1* pH, TF1* pFit) {
- pFit->SetParameter(0, 100);
- pFit->SetParameter(1, 0.001);
- pFit->SetParameter(2, 1.000);
- //pFit->SetParameter(3, 0.5);
- pH->Fit(pFit);
- };
-
- for (int iSt = 0; iSt < nSt; ++iSt) {
- FitPulls(fvph_pull_x[iSt], pPullFit);
- FitPulls(fvph_pull_y[iSt], pPullFit);
- FitPulls(fvph_pull_u[iSt], pPullFit);
- FitPulls(fvph_pull_v[iSt], pPullFit);
- FitPulls(fvph_pull_t[iSt], pPullFit);
- }
}
// ---------------------------------------------------------------------------------------------------------------------
//
-bool CbmCaInputQaSts::Check() const
+bool CbmCaInputQaSts::Check()
{
- LOG_IF(info, fVerbose > 0) << fName << ": Checking QA ...";
-
bool res = true;
int nSt = fpDetInterface->GetNtrackingStations();
@@ -107,7 +68,7 @@ bool CbmCaInputQaSts::Check() const
// ** Basic checks, available both for real and simulated data **
// **************************************************************
- // ----- Checks for mismatches in the order of stations
+ // ----- Checks for mismatches in the ordering of the stations
//
std::vector<double> vStationPos(nSt, 0.);
for (int iSt = 0; iSt < nSt; ++iSt) {
@@ -143,8 +104,8 @@ bool CbmCaInputQaSts::Check() const
res = false;
continue;
}
- int iBinMin = fvph_hit_station_delta_z[iSt]->FindBin(-cuts::kMaxDzStHitSts);
- int iBinMax = fvph_hit_station_delta_z[iSt]->FindBin(+cuts::kMaxDzStHitSts);
+ int iBinMin = fvph_hit_station_delta_z[iSt]->FindBin(-fMaxDiffZStHit);
+ int iBinMax = fvph_hit_station_delta_z[iSt]->FindBin(+fMaxDiffZStHit);
if (fvph_hit_station_delta_z[iSt]->Integral(iBinMin, iBinMax) < nHits) {
LOG_IF(error, fVerbose > 0) << fName << ": station " << iSt << " has mismatches in hit z-positions";
@@ -158,26 +119,70 @@ bool CbmCaInputQaSts::Check() const
// *******************************************************
if (IsMCUsed()) {
-
// ----- Check efficiencies
// Error is raised, if any station has integrated efficiency lower then a selected threshold.
// Integrated efficiency is estimated with fitting the efficiency curve for hit distances (r) with a uniform
// distribution in the range from kEffRangeMin to kEffRangeMax, where the efficiency is assigned to be constant
//
+ // Fit efficiency curves
LOG(info) << "-- Hit efficiency integrated over hit distance from station center";
- LOG(info) << "\tintergration range: [" << kEffRangeMin << ", " << kEffRangeMax << "] cm";
- LOG(info) << std::setw(10) << std::setfill(' ') << "St. ID" << ' ' << std::setw(12) << std::setfill(' ') << "eps.";
+ LOG(info) << "\tintegration range: [" << fEffRange[0] << ", " << fEffRange[1] << "] cm";
+
+ auto* pEffTable = MakeTable("eff_table", "Efficiency table", nSt, 2);
+ pEffTable->SetNamesOfCols({"Station ID", "Efficiency"});
+ pEffTable->SetColWidth(20);
+
for (int iSt = 0; iSt < nSt; ++iSt) {
- auto* pFitfn = (TF1*) fvpe_reco_eff_vs_r[iSt]->GetListOfFunctions()->FindObject("fitfn");
- double eff = pFitfn->GetParameter(0);
- LOG(info) << std::setw(10) << std::setfill(' ') << iSt << ' ' << std::setw(12) << std::setfill(' ') << eff;
- if (eff < cuts::kMinEff) {
- LOG(error) << fName << ": station " << iSt << " has efficiency lower threshold (" << eff << " < "
- << cuts::kMinEff << ')';
- res = false;
- }
+ auto pEff = std::unique_ptr<CbmQaEff>(fvpe_reco_eff_vs_r[iSt]->Integral(fEffRange[0], fEffRange[1]));
+ double eff = pEff->GetEfficiency(1);
+ pEffTable->SetCell(iSt, 0, iSt);
+ pEffTable->SetCell(iSt, 1, eff);
+ res = CheckRange("Hit finder efficiency in station " + std::to_string(iSt), eff, fEffThrsh, 1.000);
+ }
+
+ LOG(info) << '\n' << pEffTable->ToString(3);
+
+ // ----- Checks for residuals
+ // Check hits for potential biases from central values
+
+ // Function to fit a residual distribution, returns a structure
+ auto FitResidualsAndGetMean = [&](TH1* pHist) {
+ auto* pFit = (TF1*) gROOT->FindObject("gausn");
+ LOG_IF(fatal, !pFit) << fName << ": residuals fit function is not found";
+ double statMean = pHist->GetMean();
+ double statSigm = pHist->GetStdDev();
+ pFit->SetParameters(100., statMean, statSigm);
+ pHist->Fit(pFit, "MQ");
+ pHist->Fit(pFit, "MQ");
+ pHist->Fit(pFit, "MQ");
+ // NOTE: Several fit procedures are made to avoid empty fit results
+ std::array<double, 3> result;
+ result[0] = pFit->GetParameter(1);
+ result[1] = -pFit->GetParameter(2) * fResMeanThrsh;
+ result[2] = +pFit->GetParameter(2) * fResMeanThrsh;
+ return result;
+ };
+
+ auto* pResidualsTable = MakeTable("residuals_mean", "Residual mean values in different stations", nSt, 4);
+ pResidualsTable->SetNamesOfCols({"Station ID", "Residual(x) [cm]", "Residual(y) [cm]", "Residual(t) [ns]"});
+ pResidualsTable->SetColWidth(20);
+
+ // Fill residuals fit results and check boundaries
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ auto fitX = FitResidualsAndGetMean(fvph_res_x[iSt]);
+ auto fitY = FitResidualsAndGetMean(fvph_res_y[iSt]);
+ auto fitT = FitResidualsAndGetMean(fvph_res_t[iSt]);
+ res = CheckRange("Mean of x residual [cm] in station " + std::to_string(iSt), fitX[0], fitX[1], fitX[2]) && res;
+ res = CheckRange("Mean of y residual [cm] in station " + std::to_string(iSt), fitY[0], fitY[1], fitY[2]) && res;
+ res = CheckRange("Mean of t residual [ns] in station " + std::to_string(iSt), fitT[0], fitT[1], fitT[2]) && res;
+ pResidualsTable->SetCell(iSt, 0, iSt);
+ pResidualsTable->SetCell(iSt, 1, fitX[0]);
+ pResidualsTable->SetCell(iSt, 2, fitX[1]);
+ pResidualsTable->SetCell(iSt, 3, fitX[2]);
}
+ LOG(info) << '\n' << pResidualsTable->ToString(8);
+
// ----- Checks for pulls
// For the hit pull is defined as a ration of the difference between hit and MC-point position or time component
// values and an error of the hit position (or time) component, calculated in the same z-positions. In ideal case,
@@ -186,35 +191,45 @@ bool CbmCaInputQaSts::Check() const
// this range, QA task fails.
// TODO: Add checks for center
- // Function returns width of histogram
- auto GetPullWidth = [&](TH1* pHist) {
- if (pHist->FindObject("pullFitGausn")) {
- TF1* pFitFn = (TF1*) pHist->FindObject("pullFitGausn");
- return pFitFn->GetParameter(2);
- }
- LOG(error) << fName << ": no fit function provided for pulls fitting (" << pHist->GetName() << "), RMS returned";
- return pHist->GetRMS();
+ // Fit pull distributions
+
+ // ** Kaniadakis Gaussian distribution, gives smaller chi2 / NDF
+ //if (!gROOT->FindObject("Expk")) {
+ // new TFormula("Expk", "TMath::Power(TMath::Sqrt(1 + x[1] * x[1] * x[0] * x[0]) + x[0] * x[1], 1./x[1])");
+ //}
+ //TF1* pPullFit = new TF1("pullFitGausn", "[0] * Expk(-[2] * (x - [1]) * (x - [1]), [3])", -10., 10.);
+
+ auto FitPullsAndGetSigma = [&](TH1* pHist) {
+ auto* pFit = (TF1*) gROOT->FindObject("gausn");
+ LOG_IF(fatal, !pFit) << fName << ": pulls fit function is not found";
+ pFit->SetParameters(100, 0.001, 1.000);
+ pHist->Fit(pFit, "Q");
+ return pFit->GetParameter(2);
};
- LOG(info) << "-- Hit pull RMS:";
- LOG(info) << std::setw(10) << std::setfill(' ') << "St. ID" << ' ' << std::setw(12) << std::setfill(' ') << "RMS(x)"
- << ' ' << std::setw(12) << std::setfill(' ') << "RMS(y)" << ' ' << std::setw(12) << std::setfill(' ')
- << "RMS(u)" << ' ' << std::setw(12) << std::setfill(' ') << "RMS(v)" << ' ' << std::setw(12)
- << std::setfill(' ') << "RMS(t)";
+ auto* pPullsTable = MakeTable("pulls_rms", "Pulls std. dev. values in different stations", nSt, 4);
+ pPullsTable->SetNamesOfCols({"Station ID", "Pull(x) sigm", "Pull(y) sigm", "Pull(z) sigm"});
+ pPullsTable->SetColWidth(20);
+
for (int iSt = 0; iSt < nSt; ++iSt) {
- double rmsPullX = GetPullWidth(fvph_pull_x[iSt]);
- double rmsPullY = GetPullWidth(fvph_pull_y[iSt]);
- double rmsPullU = GetPullWidth(fvph_pull_u[iSt]);
- double rmsPullV = GetPullWidth(fvph_pull_v[iSt]);
- double rmsPullT = GetPullWidth(fvph_pull_t[iSt]);
- LOG(info) << std::setw(10) << std::setfill(' ') << iSt << ' ' << std::setw(12) << std::setfill(' ') << rmsPullX
- << ' ' << std::setw(12) << std::setfill(' ') << rmsPullY << ' ' << std::setw(12) << std::setfill(' ')
- << rmsPullU << ' ' << std::setw(12) << std::setfill(' ') << rmsPullV << ' ' << std::setw(12)
- << std::setfill(' ') << rmsPullT;
+ double rmsPullX = FitPullsAndGetSigma(fvph_pull_x[iSt]);
+ double rmsPullY = FitPullsAndGetSigma(fvph_pull_y[iSt]);
+ double rmsPullT = FitPullsAndGetSigma(fvph_pull_t[iSt]);
+
+ double rmsPullHi = 1. + fPullWidthThrsh;
+ double rmsPullLo = 1. - fPullWidthThrsh;
+
+ res = CheckRange("Rms for x pull in station " + std::to_string(iSt), rmsPullX, rmsPullLo, rmsPullHi) && res;
+ res = CheckRange("Rms for y pull in station " + std::to_string(iSt), rmsPullY, rmsPullLo, rmsPullHi) && res;
+ res = CheckRange("Rms for t pull in station " + std::to_string(iSt), rmsPullT, rmsPullLo, rmsPullHi) && res;
+ pPullsTable->SetCell(iSt, 0, iSt);
+ pPullsTable->SetCell(iSt, 1, rmsPullX);
+ pPullsTable->SetCell(iSt, 2, rmsPullY);
+ pPullsTable->SetCell(iSt, 3, rmsPullT);
}
- }
- LOG_IF(error, !res) << fName << ": QA check failed";
+ LOG(info) << '\n' << pPullsTable->ToString(3);
+ }
return res;
}
@@ -420,10 +435,10 @@ void CbmCaInputQaSts::FillHistograms()
double pMC = sqrt(pxMC * pxMC + pyMC * pyMC + pzMC * pzMC);
// MC point exit momenta
- double pxMCo = pMCPoint->GetPxOut();
- double pyMCo = pMCPoint->GetPyOut();
- double pzMCo = pMCPoint->GetPzOut();
- double pMCo = sqrt(pxMCo * pxMCo + pyMCo * pyMCo + pzMCo * pzMCo);
+ // double pxMCo = pMCPoint->GetPxOut();
+ // double pyMCo = pMCPoint->GetPyOut();
+ // double pzMCo = pMCPoint->GetPzOut();
+ // double pMCo = sqrt(pxMCo * pxMCo + pyMCo * pyMCo + pzMCo * pzMCo);
// Position and time shifted to z-coordinate of the hit
double shiftZ = zHit - zMC; // Difference between hit and point z positions
@@ -530,14 +545,13 @@ void CbmCaInputQaSts::FillHistograms()
double rMC = sqrt(xMC * xMC + yMC * yMC);
// Conditions under which point is accounted as reconstructed: point
- bool isOneHitForOnePoint = (vNofHitsPerPoint[iE][iP] == 1);
- bool isSevHitsForOnePoint = (vNofHitsPerPoint[iE][iP] > 1);
+ bool ifPointHasHits = (vNofHitsPerPoint[iE][iP] > 0);
- fvpe_reco_eff_vs_xy[iSt]->Fill(isOneHitForOnePoint, xMC, yMC);
- fvpe_reco_eff_vs_xy[nSt]->Fill(isOneHitForOnePoint, xMC, yMC);
+ fvpe_reco_eff_vs_xy[iSt]->Fill(ifPointHasHits, xMC, yMC);
+ fvpe_reco_eff_vs_xy[nSt]->Fill(ifPointHasHits, xMC, yMC);
- fvpe_reco_eff_vs_r[iSt]->Fill(isOneHitForOnePoint, rMC);
- fvpe_reco_eff_vs_r[nSt]->Fill(isOneHitForOnePoint, rMC);
+ fvpe_reco_eff_vs_r[iSt]->Fill(ifPointHasHits, rMC);
+ fvpe_reco_eff_vs_r[nSt]->Fill(ifPointHasHits, rMC);
} // loop over MC-points: end
} // loop over MC-events: end
@@ -1127,11 +1141,11 @@ InitStatus CbmCaInputQaSts::InitHistograms()
sN = (TString) "reco_eff_vs_xy" + nsuff;
sT = (TString) "Hit rec. efficiency" + tsuff + ";x_{MC} [cm];y_{MC} [cm];#epsilon";
- fvpe_reco_eff_vs_xy[iSt] = MakeEfficiency<TEfficiency>(sN, sT, 100, -50, 50, 100, -50, 50);
+ fvpe_reco_eff_vs_xy[iSt] = MakeEfficiency<CbmQaEff>(sN, sT, 100, -50, 50, 100, -50, 50);
sN = (TString) "reco_eff_vs_r" + nsuff;
sT = (TString) "Hit rec. efficiency" + tsuff + ";r_{MC} [cm];#epsilon";
- fvpe_reco_eff_vs_r[iSt] = MakeEfficiency<TEfficiency>(sN, sT, 100, -70, 70);
+ fvpe_reco_eff_vs_r[iSt] = MakeEfficiency<CbmQaEff>(sN, sT, 100, 0, 100);
}
}
return kSUCCESS;
diff --git a/reco/L1/qa/CbmCaInputQaSts.h b/reco/L1/qa/CbmCaInputQaSts.h
index 3e205c7803..343307a805 100644
--- a/reco/L1/qa/CbmCaInputQaSts.h
+++ b/reco/L1/qa/CbmCaInputQaSts.h
@@ -11,6 +11,7 @@
#ifndef CbmCaInputQaSts_h
#define CbmCaInputQaSts_h 1
+#include "CbmCaInputQaBase.h"
#include "CbmMCDataManager.h"
#include "CbmQaTask.h"
@@ -30,10 +31,10 @@ class CbmStsTrackingInterface;
class TClonesArray;
class TH1F;
class TH2F;
-class TEfficiency;
+class CbmQaEff;
/// A QA-task class, which provides assurance of MuCh hits and geometry
-class CbmCaInputQaSts : public CbmQaTask {
+class CbmCaInputQaSts : public CbmQaTask, public CbmCaInputQaBase {
public:
/// Constructor from parameters
/// \param verbose Verbose level
@@ -48,7 +49,7 @@ protected:
// ********************************************
/// Checks results of the QA and returns some flag
- bool Check() const;
+ bool Check();
/// Initializes data branches
InitStatus InitDataBranches();
@@ -62,17 +63,10 @@ protected:
/// Fills histograms per event or time-slice
void FillHistograms();
- /// Process filled histograms in the end of the run
- void AnalyzeHistograms();
-
/// De-initializes histograms
void DeInit();
private:
- // *********************
- // ** Cut variables
-
-
// *********************
// ** Private methods **
// *********************
@@ -87,8 +81,7 @@ private:
// ** Class variables **
// *********************
-
- // ** Data branches **
+ // ----- Data branches
CbmStsTrackingInterface* fpDetInterface = nullptr; ///< Instance of detector interface
CbmTimeSlice* fpTimeSlice = nullptr; ///< Pointer to current time-slice
@@ -103,11 +96,8 @@ private:
TClonesArray* fpClusterMatches = nullptr; ///< Array of hit matches
- // ** Histograms **
-
- static constexpr double kEffRangeMin = 10.; ///< Lower limit of hit distance for efficiency integration [cm]
- static constexpr double kEffRangeMax = 30.; ///< Upper limit of hit distance for efficiency integration [cm]
+ // ----- Histograms
static constexpr double kRHitX[2] = {-50., 50}; ///< Range for hit x coordinate [cm]
static constexpr double kRHitY[2] = {-50., 50}; ///< Range for hit y coordinate [cm]
static constexpr double kRHitZ[2] = {-20., 60}; ///< Range for hit z coordinate [cm]
@@ -116,7 +106,7 @@ private:
static constexpr int kNbinsZ = 800; ///< Number of bins for z coordinate axis
static constexpr double kRHitDx[2] = {-.005, .005}; ///< Range for hit x coordinate [cm]
- static constexpr double kRHitDy[2] = {-.005, .005}; ///< Range for hit y coordinate [cm]
+ static constexpr double kRHitDy[2] = {-.02, .02}; ///< Range for hit y coordinate [cm]
static constexpr double kRHitDu[2] = {-.005, .005}; ///< Range for hit u coordinate [cm]
static constexpr double kRHitDv[2] = {-.005, .005}; ///< Range for hit v coordinate [cm]
static constexpr double kRHitDt[2] = {-10., 10.}; ///< Range for hit time [ns]
@@ -133,18 +123,16 @@ private:
static constexpr double kRPullV[2] = {-10., 10.}; ///< Range for pull of y coordinate
static constexpr double kRPullT[2] = {-5., 5.}; ///< Range for pull of time
-
// NOTE: the last element of each vector stands for integral distribution over all stations
- // hit occupancy
+ // Hit occupancy
std::vector<TH2F*> fvph_hit_ypos_vs_xpos; ///< hit ypos vs xpos in different stations
std::vector<TH2F*> fvph_hit_xpos_vs_zpos; ///< hit xpos vs zpos in different stations
std::vector<TH2F*> fvph_hit_ypos_vs_zpos; ///< hit ypos vs zpos in different stations
- // difference between hit and station z positions
std::vector<TH1F*> fvph_hit_station_delta_z; ///< Difference between station and hit z positions [cm]
- // hit errors
+ // Hit errors
std::vector<TH1F*> fvph_hit_dx;
std::vector<TH1F*> fvph_hit_dy;
std::vector<TH1F*> fvph_hit_du;
@@ -154,7 +142,6 @@ private:
// MC points occupancy
std::vector<TH1F*> fvph_n_points_per_hit; ///< number of points per one hit in different stations
-
std::vector<TH2F*> fvph_point_ypos_vs_xpos; ///< point ypos vs xpos in different stations
std::vector<TH2F*> fvph_point_xpos_vs_zpos; ///< point xpos vs zpos in different stations
std::vector<TH2F*> fvph_point_ypos_vs_zpos; ///< point ypos vs zpos in different stations
@@ -188,14 +175,8 @@ private:
std::vector<TH2F*> fvph_pull_t_vs_t; ///< pull for t coordinate in different stations
// Hit efficiencies
- std::vector<TEfficiency*> fvpe_reco_eff_vs_xy; ///< Efficiency of hit reconstruction vs. x and y coordinates of MC
- std::vector<TEfficiency*> fvpe_reco_eff_vs_r; ///< Efficiency of hit reconstruction vs. distance from center
-
- /// Kaniadakis exponent
- static Double_t Expk(Double_t x, Double_t k)
- {
- return TMath::Power((TMath::Sqrt(1 + k * k * x * x) + k * x), 1. / k);
- }
+ std::vector<CbmQaEff*> fvpe_reco_eff_vs_xy; ///< Efficiency of hit reconstruction vs. x and y coordinates of MC
+ std::vector<CbmQaEff*> fvpe_reco_eff_vs_r; ///< Efficiency of hit reconstruction vs. distance from center
};
#endif // CbmCaInputQaSts_h
diff --git a/reco/L1/qa/CbmCaInputQaTof.cxx b/reco/L1/qa/CbmCaInputQaTof.cxx
new file mode 100644
index 0000000000..cb0e89ae7f
--- /dev/null
+++ b/reco/L1/qa/CbmCaInputQaTof.cxx
@@ -0,0 +1,1108 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// \file CbmCaInputQaTof.cxx
+/// \date 30.01.2023
+/// \brief QA-task for CA tracking input from TOF detector (implementation)
+/// \author S.Zharko <s.zharko@gsi.de>
+
+#include "CbmCaInputQaTof.h"
+
+#include "CbmAddress.h"
+#include "CbmLink.h"
+#include "CbmMCDataArray.h"
+#include "CbmMCEventList.h"
+#include "CbmMCTrack.h"
+#include "CbmMatch.h"
+#include "CbmQaCanvas.h"
+#include "CbmQaEff.h"
+#include "CbmTimeSlice.h"
+#include "CbmTofHit.h"
+#include "CbmTofInteraction.h"
+#include "CbmTofPoint.h"
+#include "CbmTofTrackingInterface.h"
+
+#include "FairRootManager.h"
+#include "Logger.h"
+
+#include "TBox.h"
+#include "TClonesArray.h"
+#include "TEllipse.h"
+#include "TF1.h"
+#include "TFormula.h"
+#include "TGraphAsymmErrors.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TMath.h"
+#include "TStyle.h"
+
+#include <algorithm>
+#include <fstream>
+#include <iomanip>
+#include <numeric>
+#include <set> // TMP!!!!
+
+#include "CaToolsLinkKey.h"
+#include "L1Constants.h"
+
+ClassImp(CbmCaInputQaTof);
+
+namespace phys = L1Constants::phys; // from L1Constants.h
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+CbmCaInputQaTof::CbmCaInputQaTof(int verbose, bool isMCUsed) : CbmQaTask("CbmCaInputQaTof", "catof", verbose, isMCUsed)
+{
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+bool CbmCaInputQaTof::Check()
+{
+ bool res = true;
+
+ int nSt = fpDetInterface->GetNtrackingStations();
+
+
+ // **************************************************************
+ // ** Basic checks, available both for real and simulated data **
+ // **************************************************************
+
+ // ----- Checks for mismatches in the ordering of the stations
+ //
+ std::vector<double> vStationPos(nSt, 0.);
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ vStationPos[iSt] = fpDetInterface->GetZ(iSt);
+ }
+
+ if (!std::is_sorted(vStationPos.cbegin(), vStationPos.cend(), [](int l, int r) { return l <= r; })) {
+ if (fVerbose > 0) {
+ LOG(error) << fName << ": stations are ordered improperly along the beam axis:";
+ for (auto zPos : vStationPos) {
+ LOG(error) << "\t- " << zPos;
+ }
+ }
+ res = false;
+ }
+
+ // ----- Checks for mismatch between station and hit z positions
+ // The purpose of this block is to be ensured, that hits belong to the correct tracking station. For each tracking
+ // station a unified position z-coordinate is defined, which generally differs from the corresponding positions of
+ // reconstructed hits. This is due to non-regular position along the beam axis for each detector sensor. Thus, the
+ // positions of the hits along z-axis are distributed relatively to the position of the station in some range.
+ // If hits goes out this range, it can signal about a mismatch between hits and geometry. For limits of the range
+ // one should select large enough values to cover the difference described above and in the same time small enough
+ // to avoid overlaps with the neighboring stations.
+ // For each station, a distribution of z_{hit} - z_{st} is integrated over the defined range and scaled by the
+ // total number of entries to the distribution. If this value is smaller then unity, then some of the hits belong to
+ // another station.
+ //
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ int nHits = fvph_hit_station_delta_z[iSt]->GetEntries();
+ if (!nHits) {
+ LOG_IF(error, fVerbose > 0) << fName << ": station " << iSt << " does not have hits";
+ res = false;
+ continue;
+ }
+ int iBinMin = fvph_hit_station_delta_z[iSt]->FindBin(-fMaxDiffZStHit);
+ int iBinMax = fvph_hit_station_delta_z[iSt]->FindBin(+fMaxDiffZStHit);
+
+ if (fvph_hit_station_delta_z[iSt]->Integral(iBinMin, iBinMax) < nHits) {
+ LOG_IF(error, fVerbose > 0) << fName << ": station " << iSt << " has mismatches in hit z-positions";
+ res = false;
+ }
+ }
+
+
+ // *******************************************************
+ // ** Additional checks, if MC information is available **
+ // *******************************************************
+
+ if (IsMCUsed()) {
+ // ----- Check efficiencies
+ // Error is raised, if any station has integrated efficiency lower then a selected threshold.
+ // Integrated efficiency is estimated with fitting the efficiency curve for hit distances (r) with a uniform
+ // distribution in the range from kEffRangeMin to kEffRangeMax, where the efficiency is assigned to be constant
+ //
+ // Fit efficiency curves
+ LOG(info) << "-- Hit efficiency integrated over hit distance from station center";
+ LOG(info) << "\tintegration range: [" << fEffRange[0] << ", " << fEffRange[1] << "] cm";
+
+ auto* pEffTable = MakeTable("eff_table", "Efficiency table", nSt, 2);
+ pEffTable->SetNamesOfCols({"Station ID", "Efficiency"});
+ pEffTable->SetColWidth(20);
+
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ auto pEff = std::unique_ptr<CbmQaEff>(fvpe_reco_eff_vs_r[iSt]->Integral(fEffRange[0], fEffRange[1]));
+ double eff = pEff->GetEfficiency(1);
+ pEffTable->SetCell(iSt, 0, iSt);
+ pEffTable->SetCell(iSt, 1, eff);
+ res = CheckRange("Hit finder efficiency in station " + std::to_string(iSt), eff, fEffThrsh, 1.000);
+ }
+
+ LOG(info) << '\n' << pEffTable->ToString(3);
+
+ // ----- Checks for residuals
+ // Check hits for potential biases from central values
+
+ // Function to fit a residual distribution, returns a structure
+ auto FitResidualsAndGetMean = [&](TH1* pHist) {
+ auto* pFit = (TF1*) gROOT->FindObject("gausn");
+ LOG_IF(fatal, !pFit) << fName << ": residuals fit function is not found";
+ double statMean = pHist->GetMean();
+ double statSigm = pHist->GetStdDev();
+ pFit->SetParameters(100., statMean, statSigm);
+ pHist->Fit(pFit, "MQ");
+ pHist->Fit(pFit, "MQ");
+ pHist->Fit(pFit, "MQ");
+ // NOTE: Several fit procedures are made to avoid empty fit results
+ std::array<double, 3> result;
+ result[0] = pFit->GetParameter(1);
+ result[1] = -pFit->GetParameter(2) * fResMeanThrsh;
+ result[2] = +pFit->GetParameter(2) * fResMeanThrsh;
+ return result;
+ };
+
+ auto* pResidualsTable = MakeTable("residuals_mean", "Residual mean values in different stations", nSt, 4);
+ pResidualsTable->SetNamesOfCols({"Station ID", "Residual(x) [cm]", "Residual(y) [cm]", "Residual(t) [ns]"});
+ pResidualsTable->SetColWidth(20);
+
+ // Fill residuals fit results and check boundaries
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ auto fitX = FitResidualsAndGetMean(fvph_res_x[iSt]);
+ auto fitY = FitResidualsAndGetMean(fvph_res_y[iSt]);
+ auto fitT = FitResidualsAndGetMean(fvph_res_t[iSt]);
+ res = CheckRange("Mean of x residual [cm] in station " + std::to_string(iSt), fitX[0], fitX[1], fitX[2]) && res;
+ res = CheckRange("Mean of y residual [cm] in station " + std::to_string(iSt), fitY[0], fitY[1], fitY[2]) && res;
+ res = CheckRange("Mean of t residual [ns] in station " + std::to_string(iSt), fitT[0], fitT[1], fitT[2]) && res;
+ pResidualsTable->SetCell(iSt, 0, iSt);
+ pResidualsTable->SetCell(iSt, 1, fitX[0]);
+ pResidualsTable->SetCell(iSt, 2, fitX[1]);
+ pResidualsTable->SetCell(iSt, 3, fitX[2]);
+ }
+
+ LOG(info) << '\n' << pResidualsTable->ToString(8);
+
+ // ----- Checks for pulls
+ // For the hit pull is defined as a ration of the difference between hit and MC-point position or time component
+ // values and an error of the hit position (or time) component, calculated in the same z-positions. In ideal case,
+ // when the resolutions are well determined for detecting stations, the pull distributions should have a RMS equal
+ // to unity. Here we allow a RMS of the pull distributions to be varied in a predefined range. If the RMS runs out
+ // this range, QA task fails.
+ // TODO: Add checks for center
+
+ // Fit pull distributions
+
+ // ** Kaniadakis Gaussian distribution, gives smaller chi2 / NDF
+ //if (!gROOT->FindObject("Expk")) {
+ // new TFormula("Expk", "TMath::Power(TMath::Sqrt(1 + x[1] * x[1] * x[0] * x[0]) + x[0] * x[1], 1./x[1])");
+ //}
+ //TF1* pPullFit = new TF1("pullFitGausn", "[0] * Expk(-[2] * (x - [1]) * (x - [1]), [3])", -10., 10.);
+
+ auto FitPullsAndGetSigma = [&](TH1* pHist) {
+ auto* pFit = (TF1*) gROOT->FindObject("gausn");
+ LOG_IF(fatal, !pFit) << fName << ": pulls fit function is not found";
+ pFit->SetParameters(100, 0.001, 1.000);
+ pHist->Fit(pFit, "Q");
+ return pFit->GetParameter(2);
+ };
+
+ auto* pPullsTable = MakeTable("pulls_rms", "Pulls std. dev. values in different stations", nSt, 4);
+ pPullsTable->SetNamesOfCols({"Station ID", "Pull(x) sigm", "Pull(y) sigm", "Pull(z) sigm"});
+ pPullsTable->SetColWidth(20);
+
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ double rmsPullX = FitPullsAndGetSigma(fvph_pull_x[iSt]);
+ double rmsPullY = FitPullsAndGetSigma(fvph_pull_y[iSt]);
+ double rmsPullT = FitPullsAndGetSigma(fvph_pull_t[iSt]);
+
+ double rmsPullHi = 1. + fPullWidthThrsh;
+ double rmsPullLo = 1. - fPullWidthThrsh;
+
+ res = CheckRange("Rms for x pull in station " + std::to_string(iSt), rmsPullX, rmsPullLo, rmsPullHi) && res;
+ res = CheckRange("Rms for y pull in station " + std::to_string(iSt), rmsPullY, rmsPullLo, rmsPullHi) && res;
+ res = CheckRange("Rms for t pull in station " + std::to_string(iSt), rmsPullT, rmsPullLo, rmsPullHi) && res;
+ pPullsTable->SetCell(iSt, 0, iSt);
+ pPullsTable->SetCell(iSt, 1, rmsPullX);
+ pPullsTable->SetCell(iSt, 2, rmsPullY);
+ pPullsTable->SetCell(iSt, 3, rmsPullT);
+ }
+
+ LOG(info) << '\n' << pPullsTable->ToString(3);
+ }
+
+ return res;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void CbmCaInputQaTof::DeInit()
+{
+ // Vectors with pointers to histograms
+ fvph_hit_ypos_vs_xpos.clear();
+ fvph_hit_xpos_vs_zpos.clear();
+ fvph_hit_ypos_vs_zpos.clear();
+
+ fvph_hit_station_delta_z.clear();
+
+ fvph_hit_dx.clear();
+ fvph_hit_dy.clear();
+ fvph_hit_dt.clear();
+
+ fvph_n_points_per_hit.clear();
+
+ fvph_point_ypos_vs_xpos.clear();
+ fvph_point_xpos_vs_zpos.clear();
+ fvph_point_ypos_vs_zpos.clear();
+
+ fvph_point_hit_delta_z.clear();
+
+ fvph_res_x.clear();
+ fvph_res_y.clear();
+ fvph_res_t.clear();
+
+ fvph_pull_x.clear();
+ fvph_pull_y.clear();
+ fvph_pull_t.clear();
+
+ fvph_res_x_vs_x.clear();
+ fvph_res_y_vs_y.clear();
+ fvph_res_t_vs_t.clear();
+
+ fvph_pull_x_vs_x.clear();
+ fvph_pull_y_vs_y.clear();
+ fvph_pull_t_vs_t.clear();
+
+ fvpe_reco_eff_vs_xy.clear();
+ fvpe_reco_eff_vs_r.clear();
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void CbmCaInputQaTof::FillHistograms()
+{
+ int nSt = fpDetInterface->GetNtrackingStations();
+ int nHits = fpHits->GetEntriesFast();
+
+ // ----- Fill TOF interactions
+ std::shared_ptr<TClonesArray> pMCInteractions = nullptr; // Array of MC interactions
+ std::shared_ptr<TClonesArray> pHitInterMatches = nullptr; // Array of hit matches to MC interactions
+ if (IsMCUsed()) {
+ pMCInteractions = std::make_shared<TClonesArray>("CbmTofInteraction");
+ pHitInterMatches = std::make_shared<TClonesArray>("CbmMatch", nHits);
+ FillInteractions(pMCInteractions, pHitInterMatches);
+ } // IsMCUsed
+
+ std::vector<int> vNofHitsPerInteraction; // Number of hits per MC point in TS
+
+ if (IsMCUsed()) { vNofHitsPerInteraction.resize(pMCInteractions->GetEntriesFast()); }
+
+ for (int iH = 0; iH < nHits; ++iH) {
+ const auto* pHit = dynamic_cast<const CbmTofHit*>(fpHits->At(iH));
+ LOG_IF(fatal, !pHit) << fName << ": hit with iH = " << iH << " is not an CbmTofHit (dynamic cast failed)";
+
+ // FIXME: Is this cut still needed??
+ int32_t address = pHit->GetAddress();
+ if (0x00202806 == address || 0x00002806 == address) { continue; }
+
+ // *************************
+ // ** Reconstructed hit QA
+
+ // Hit station geometry info
+ int iSt = fpDetInterface->GetTrackingStationIndex(pHit);
+ LOG_IF(fatal, iSt < 0 || iSt >= nSt) << fName << ": index of station (" << iSt << ") is out of range "
+ << "for hit with id = " << iH;
+
+ // Hit position
+ double xHit = pHit->GetX();
+ double yHit = pHit->GetY();
+ double zHit = pHit->GetZ();
+
+ double dxHit = pHit->GetDx();
+ double dyHit = pHit->GetDy();
+ //double dxyHit = pHit->GetDxy();
+
+ // Hit time
+ double tHit = pHit->GetTime();
+ double dtHit = pHit->GetTimeError();
+
+ fvph_hit_ypos_vs_xpos[iSt]->Fill(xHit, yHit);
+ fvph_hit_xpos_vs_zpos[iSt]->Fill(zHit, xHit);
+ fvph_hit_ypos_vs_zpos[iSt]->Fill(zHit, yHit);
+
+ fvph_hit_station_delta_z[iSt]->Fill(zHit - fpDetInterface->GetZ(iSt));
+
+ fvph_hit_dx[iSt]->Fill(dxHit);
+ fvph_hit_dy[iSt]->Fill(dyHit);
+ fvph_hit_dt[iSt]->Fill(dtHit);
+
+ fvph_hit_ypos_vs_xpos[nSt]->Fill(xHit, yHit);
+ fvph_hit_xpos_vs_zpos[nSt]->Fill(zHit, xHit);
+ fvph_hit_ypos_vs_zpos[nSt]->Fill(zHit, yHit);
+ fvph_hit_dx[nSt]->Fill(dxHit);
+ fvph_hit_dy[nSt]->Fill(dyHit);
+ fvph_hit_dt[nSt]->Fill(dtHit);
+
+
+ // **********************
+ // ** MC information QA
+
+ if (IsMCUsed()) {
+ CbmMatch* pHitMatch = dynamic_cast<CbmMatch*>(pHitInterMatches->At(iH));
+ LOG_IF(fatal, !pHitMatch) << fName << ": match object not found for hit ID " << iH;
+
+ // NOTE: Here we treat interactions simply as MC points
+ // Evaluate number of hits per one MC point
+ int nMCpoints = 0; // Number of MC points for this hit
+ int nLinks = pHitMatch->GetNofLinks();
+ for (int iLink = 0; iLink < nLinks; ++iLink) {
+ const auto& link = pHitMatch->GetLink(iLink);
+
+ int iP = link.GetIndex(); // Index of MC interaction
+
+ // Skip noisy links
+ if (iP < 0) { continue; }
+
+ ++nMCpoints;
+
+ LOG_IF(fatal, iP >= (int) vNofHitsPerInteraction.size())
+ << fName << ": index of MC point is out of range (hit id = " << iH << ", link id = " << iLink
+ << ", point id = " << iP << ')';
+ vNofHitsPerInteraction[iP]++;
+ }
+
+ fvph_n_points_per_hit[iSt]->Fill(nMCpoints);
+ fvph_n_points_per_hit[nSt]->Fill(nMCpoints);
+
+ if (nMCpoints != 1) { continue; } // ??
+
+ // The best link to in the match (probably, the cut on nMCpoints is meaningless)
+ const auto& bestPointLink = pHitMatch->GetMatchedLink();
+
+ // Skip noisy links
+ if (bestPointLink.GetIndex() < 0) { continue; }
+
+ // Point matched by the best link
+ const auto* pMCPoint = dynamic_cast<const CbmTofInteraction*>(pMCInteractions->At(bestPointLink.GetIndex()));
+ LOG_IF(fatal, !pMCPoint) << fName << ": MC point object does not exist for hit " << iH;
+
+ // MC track for this point
+ CbmLink bestTrackLink = bestPointLink;
+ bestTrackLink.SetIndex(pMCPoint->GetTrackID());
+ const auto* pMCTrack = dynamic_cast<const CbmMCTrack*>(fpMCTracks->Get(bestTrackLink));
+ LOG_IF(fatal, !pMCTrack) << fName << ": MC track object does not exist for hit " << iH << " and link: ";
+
+ double t0MC = fpMCEventList->GetEventTime(bestPointLink);
+ LOG_IF(fatal, t0MC < 0) << fName << ": MC time zero is lower then 0 ns: " << t0MC;
+
+ // ----- MC point properties
+ //
+ double mass = pMCTrack->GetMass();
+ //double charge = pMCTrack->GetCharge();
+ //double pdg = pMCTrack->GetPdgCode();
+
+ // Entrance position and time
+ double xMC = pMCPoint->GetX();
+ double yMC = pMCPoint->GetY();
+ double zMC = pMCPoint->GetZ();
+ double tMC = pMCPoint->GetTime() + t0MC;
+
+ // MC point entrance momenta
+ double pxMC = pMCPoint->GetPx();
+ double pyMC = pMCPoint->GetPy();
+ double pzMC = pMCPoint->GetPz();
+ double pMC = sqrt(pxMC * pxMC + pyMC * pyMC + pzMC * pzMC);
+
+ // Position and time shifted to z-coordinate of the hit
+ double shiftZ = zHit - zMC; // Difference between hit and point z positions
+ double xMCs = xMC + shiftZ * pxMC / pzMC;
+ double yMCs = yMC + shiftZ * pyMC / pzMC;
+ double tMCs = tMC + shiftZ / (pzMC * phys::kSpeedOfLight) * sqrt(mass * mass + pMC * pMC);
+
+ // Residuals
+ double xRes = xHit - xMCs;
+ double yRes = yHit - yMCs;
+ double tRes = tHit - tMCs;
+
+ // ------ Cuts
+
+ if (std::fabs(pMCPoint->GetPz()) < fMinMomentum) { continue; } // CUT ON MINIMUM MOMENTUM
+ //if (pMCo < cuts::kMinP) { continue; } // Momentum threshold
+
+
+ fvph_point_ypos_vs_xpos[iSt]->Fill(xMC, yMC);
+ fvph_point_xpos_vs_zpos[iSt]->Fill(zMC, xMC);
+ fvph_point_ypos_vs_zpos[iSt]->Fill(zMC, yMC);
+
+ fvph_point_hit_delta_z[iSt]->Fill(shiftZ);
+
+ fvph_res_x[iSt]->Fill(xRes);
+ fvph_res_y[iSt]->Fill(yRes);
+ fvph_res_t[iSt]->Fill(tRes);
+
+ fvph_pull_x[iSt]->Fill(xRes / dxHit);
+ fvph_pull_y[iSt]->Fill(yRes / dyHit);
+ fvph_pull_t[iSt]->Fill(tRes / dtHit);
+
+ fvph_res_x_vs_x[iSt]->Fill(xHit, xRes);
+ fvph_res_y_vs_y[iSt]->Fill(yHit, yRes);
+ fvph_res_t_vs_t[iSt]->Fill(tHit, tRes);
+
+ fvph_pull_x_vs_x[iSt]->Fill(xHit, xRes / dxHit);
+ fvph_pull_y_vs_y[iSt]->Fill(yHit, yRes / dyHit);
+ fvph_pull_t_vs_t[iSt]->Fill(tHit, tRes / dtHit);
+
+ fvph_point_ypos_vs_xpos[nSt]->Fill(xMC, yMC);
+ fvph_point_xpos_vs_zpos[nSt]->Fill(zMC, xMC);
+ fvph_point_ypos_vs_zpos[nSt]->Fill(zMC, yMC);
+
+ fvph_point_hit_delta_z[nSt]->Fill(shiftZ);
+
+ fvph_res_x[nSt]->Fill(xRes);
+ fvph_res_y[nSt]->Fill(yRes);
+ fvph_res_t[nSt]->Fill(tRes);
+
+ fvph_pull_x[nSt]->Fill(xRes / dxHit);
+ fvph_pull_y[nSt]->Fill(yRes / dyHit);
+ fvph_pull_t[nSt]->Fill(tRes / dtHit);
+
+ fvph_res_x_vs_x[nSt]->Fill(xHit, xRes);
+ fvph_res_y_vs_y[nSt]->Fill(yHit, yRes);
+ fvph_res_t_vs_t[nSt]->Fill(tHit, tRes);
+
+ fvph_pull_x_vs_x[nSt]->Fill(xHit, xRes / dxHit);
+ fvph_pull_y_vs_y[nSt]->Fill(yHit, yRes / dyHit);
+ fvph_pull_t_vs_t[nSt]->Fill(tHit, tRes / dtHit);
+ }
+ } // loop over hits: end
+
+ // Fill hit reconstruction efficiencies
+ if (IsMCUsed()) {
+ int nPoints = pMCInteractions->GetEntriesFast();
+
+ for (int iP = 0; iP < nPoints; ++iP) {
+ const auto* pMCPoint = dynamic_cast<const CbmTofInteraction*>(pMCInteractions->At(iP));
+ LOG_IF(fatal, !pMCPoint) << fName << ": MC point does not exist for interaction " << iP;
+ //int address = pMCPoint->GetDetectorID();
+ int iSt = GetStationID(pMCPoint);
+ LOG_IF(fatal, iSt < 0 || iSt >= nSt)
+ << fName << ": MC point for index = " << iP << " has wrong station index: iSt = " << iSt;
+
+ double xMC = pMCPoint->GetX();
+ double yMC = pMCPoint->GetY();
+ double rMC = sqrt(xMC * xMC + yMC * yMC);
+
+ // Conditions under which point is accounted as reconstructed: point
+ bool ifPointHasHits = (vNofHitsPerInteraction[iP] > 0);
+
+ fvpe_reco_eff_vs_xy[iSt]->Fill(ifPointHasHits, xMC, yMC);
+ fvpe_reco_eff_vs_xy[nSt]->Fill(ifPointHasHits, xMC, yMC);
+
+ fvpe_reco_eff_vs_r[iSt]->Fill(ifPointHasHits, rMC);
+ fvpe_reco_eff_vs_r[nSt]->Fill(ifPointHasHits, rMC);
+ } // iP
+ } // MC is used: end
+
+ // Clear TClonesArray objects
+ if (IsMCUsed()) {
+ pHitInterMatches->Clear();
+ pMCInteractions->Clear();
+ } // IsMCUsed
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void CbmCaInputQaTof::FillInteractions(std::shared_ptr<TClonesArray>& pvInters,
+ std::shared_ptr<TClonesArray>& pvMatches)
+{
+ std::unordered_map<ca::tools::LinkKey, int> mPointToInter; // Map of point (defined by a link key) to interaction
+
+ int nHits = fpHits->GetEntriesFast();
+ int nMCevents = (IsMCUsed()) ? fpMCEventList->GetNofEvents() : -1;
+
+ // ----- Fill array of interactions
+ int iInter = -1; // Index of interaction (global over all MC events)
+ for (int iE = 0; iE < nMCevents; ++iE) {
+ int iFile = fpMCEventList->GetFileIdByIndex(iE);
+ int iEvent = fpMCEventList->GetEventIdByIndex(iE);
+ int nMCPoints = fpMCPoints->Size(iFile, iEvent);
+
+ int iDprev = -1; // detector ID for previous point
+ int iTprev = -1; // track ID for previous point
+ for (int iP = 0; iP < nMCPoints; ++iP) {
+ const auto* pMCPoint = dynamic_cast<const CbmTofPoint*>(fpMCPoints->Get(iFile, iEvent, iP));
+ LOG_IF(fatal, !pMCPoint) << fName << ": MC point is not defined for link (f,e,i) = " << iFile << ' ' << iEvent
+ << ' ' << iP;
+
+ int iDcurr = pMCPoint->GetDetectorID(); // Current detector ID
+ int iTcurr = pMCPoint->GetTrackID(); // Current track ID
+ if (iDcurr != iDprev || iTcurr != iTprev) {
+ iInter++;
+ iDprev = iDcurr;
+ iTprev = iTcurr;
+ new ((*pvInters)[iInter]) CbmTofInteraction(); // Add empty interaction object
+ }
+
+ // Update current intercation with point, and updtate the map
+ auto* pCurrInter = static_cast<CbmTofInteraction*>(pvInters->At(iInter));
+ pCurrInter->AddPoint(pMCPoint);
+ mPointToInter[ca::tools::LinkKey(iP, iEvent, iFile)] = iInter;
+ } // iP
+ } // iE
+
+ // ----- Match interactions with hits
+ for (int iH = 0; iH < nHits; ++iH) {
+ const auto* pHit = dynamic_cast<const CbmTofHit*>(fpHits->At(iH));
+ LOG_IF(fatal, !pHit) << fName << ": Hit is not defined for iH = " << iH;
+
+ // FIXME: Is this cut still needed??
+ int32_t address = pHit->GetAddress();
+ if (0x00202806 == address || 0x00002806 == address) { continue; }
+
+ auto* pHitInterMatch = new ((*pvMatches)[iH]) CbmMatch();
+
+ // Collect create interaction links from point links and fill the match
+ auto* pHitPointMatch = dynamic_cast<CbmMatch*>(fpHitMatches->At(iH));
+ int nPointLinks = pHitPointMatch->GetNofLinks();
+ for (int iL = 0; iL < nPointLinks; ++iL) {
+ const auto& link = pHitPointMatch->GetLink(iL);
+ int iFile = link.GetFile();
+ int iEvent = link.GetEntry();
+ int iP = link.GetIndex();
+ double weight = link.GetWeight();
+ int iInteraction = mPointToInter[ca::tools::LinkKey(iP, iEvent, iFile)];
+
+ // Create a TOF interaction link with a proper index
+ pHitInterMatch->AddLink(weight, iInteraction, iEvent, iFile);
+
+ // Update interaction with information from point matched to the hit
+ // Since the averaged interaction in general may not be on the MC-trajectory of the particle,
+ // we decided to shift the interaction to the matched point. Here we assume, that there is only one
+ // matched point from the interaction.
+ auto* pPoint = static_cast<CbmTofPoint*>(fpMCPoints->Get(iFile, iEvent, iP));
+ auto* pInter = static_cast<CbmTofInteraction*>(pvInters->At(iInteraction));
+ pInter->SetFromPoint(pPoint);
+ }
+ } // iH
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+int CbmCaInputQaTof::GetStationID(const CbmTofPoint* pPoint) const
+{
+ double dist = 1000.; // NOTE: arbitrary large number
+ int iStSelect = -1;
+ // We select the station, which center is closest to the MC point
+ for (int iSt = 0; iSt < fpDetInterface->GetNtrackingStations(); ++iSt) {
+ if (std::fabs(pPoint->GetZ() - fpDetInterface->GetZ(iSt)) < dist) {
+ dist = std::fabs(pPoint->GetZ() - fpDetInterface->GetZ(iSt));
+ iStSelect = iSt;
+ }
+ }
+ return iStSelect;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+InitStatus CbmCaInputQaTof::InitDataBranches()
+{
+ // STS tracking detector interface
+ fpDetInterface = CbmTofTrackingInterface::Instance();
+
+ LOG_IF(fatal, !fpDetInterface) << "\033[1;31m" << fName << ": tracking detector interface is undefined\033[0m";
+
+ // FairRootManager
+ auto* pFairRootManager = FairRootManager::Instance();
+ LOG_IF(fatal, !pFairRootManager) << "\033[1;31m" << fName << ": FairRootManager instance is a null pointer\033[0m";
+
+ // Time-slice
+ fpTimeSlice = dynamic_cast<CbmTimeSlice*>(pFairRootManager->GetObject("TimeSlice."));
+ LOG_IF(fatal, !fpTimeSlice) << "\033[1;31m" << fName << ": time-slice branch is not found\033[0m";
+
+ // ----- Reconstructed data-branches initialization
+
+ // Hits container
+ fpHits = dynamic_cast<TClonesArray*>(pFairRootManager->GetObject("TofHit"));
+ LOG_IF(fatal, !fpHits) << "\033[1;31m" << fName << ": container of reconstructed hits in TOF is not found\033[0m";
+
+ // ----- MC-information branches initialization
+ if (IsMCUsed()) {
+ // MC manager
+ fpMCDataManager = dynamic_cast<CbmMCDataManager*>(pFairRootManager->GetObject("MCDataManager"));
+ LOG_IF(fatal, !fpMCDataManager) << "\033[1;31m" << fName << ": MC data manager branch is not found\033[0m";
+
+ // MC event list
+ fpMCEventList = dynamic_cast<CbmMCEventList*>(pFairRootManager->GetObject("MCEventList."));
+ LOG_IF(fatal, !fpMCEventList) << "\033[1;31m" << fName << ": MC event list branch is not found\033[0m";
+
+ // MC tracks
+ fpMCTracks = fpMCDataManager->InitBranch("MCTrack");
+ LOG_IF(fatal, !fpMCTracks) << "\033[1;31m" << fName << ": MC track branch is not found\033[0m";
+
+ // MC points
+ fpMCPoints = fpMCDataManager->InitBranch("TofPoint");
+ LOG_IF(fatal, !fpMCTracks) << "\033[1;31m" << fName << ": MC point branch is not found for TOF\033[0m";
+
+ // Hit matches
+ fpHitMatches = dynamic_cast<TClonesArray*>(pFairRootManager->GetObject("TofHitMatch"));
+ LOG_IF(fatal, !fpHitMatches) << "\033[1;31m]" << fName << ": hit match branch is not found for TOF\033[0m";
+ }
+
+ return kSUCCESS;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+InitStatus CbmCaInputQaTof::InitCanvases()
+{
+ gStyle->SetOptFit(1);
+ int nSt = fpDetInterface->GetNtrackingStations();
+
+ // ********************
+ // ** Drawing options
+
+ // Contours
+ constexpr auto contColor = kOrange + 7;
+ constexpr auto contWidth = 2; // Line width
+ constexpr auto contStyle = 2; // Line style
+ constexpr auto contFill = 0; // Fill style
+
+ // *********************************
+ // ** Hit occupancies
+
+ // ** Occupancy in XY plane
+ auto* pc_hit_ypos_vs_xpos = MakeCanvas<CbmQaCanvas>("hit_ypos_vs_xpos", "", 1600, 800);
+ pc_hit_ypos_vs_xpos->Divide2D(nSt);
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ pc_hit_ypos_vs_xpos->cd(iSt + 1);
+ fvph_hit_ypos_vs_xpos[iSt]->DrawCopy("colz", "");
+
+
+ // Reference size of the station saved to the detector interface
+ auto* pInnerCircle = new TEllipse(0., 0., fpDetInterface->GetRmin(iSt));
+ pInnerCircle->SetLineWidth(contWidth);
+ pInnerCircle->SetLineStyle(contStyle);
+ pInnerCircle->SetLineColor(contColor);
+ pInnerCircle->SetFillStyle(contFill);
+ pInnerCircle->Draw("SAME");
+
+ auto* pOuterCircle = new TEllipse(0., 0., fpDetInterface->GetRmax(iSt));
+ pOuterCircle->SetLineWidth(contWidth);
+ pOuterCircle->SetLineStyle(contStyle);
+ pOuterCircle->SetLineColor(contColor);
+ pOuterCircle->SetFillStyle(contFill);
+ pOuterCircle->Draw("SAME");
+
+ // Square boarders of station, which are used for the field approximation
+ double stXmin = +fpDetInterface->GetXmax(iSt);
+ double stXmax = -fpDetInterface->GetXmax(iSt);
+ double stYmin = -fpDetInterface->GetYmax(iSt);
+ double stYmax = +fpDetInterface->GetYmax(iSt);
+
+ auto* pBox = new TBox(stXmin, stYmin, stXmax, stYmax);
+ pBox->SetLineWidth(contWidth);
+ pBox->SetLineStyle(contStyle);
+ pBox->SetLineColor(contColor);
+ pBox->SetFillStyle(contFill);
+ pBox->Draw("SAME");
+ }
+
+ // ----- Occupancy projections
+ auto* pc_hit_xpos = MakeCanvas<CbmQaCanvas>("hit_xpos", "", 1400, 800);
+ pc_hit_xpos->Divide2D(nSt);
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ pc_hit_xpos->cd(iSt + 1);
+ auto* phProj = fvph_hit_ypos_vs_xpos[iSt]->ProjectionY();
+ phProj->DrawCopy();
+ }
+
+ auto* pc_hit_ypos = MakeCanvas<CbmQaCanvas>("hit_ypos", "", 1400, 800);
+ pc_hit_ypos->Divide2D(nSt);
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ pc_hit_ypos->cd(iSt + 1);
+ auto* phProj = fvph_hit_ypos_vs_xpos[iSt]->ProjectionX();
+ phProj->DrawCopy();
+ }
+
+
+ // ** Occupancy in XZ plane
+ MakeCanvas<CbmQaCanvas>("hit_xpos_vs_zpos", "", 600, 400);
+ fvph_hit_xpos_vs_zpos[nSt]->DrawCopy("colz", "");
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ // Station positions in detector IFS
+ double stZmin = fpDetInterface->GetZ(iSt) - 0.5 * fpDetInterface->GetThickness(iSt);
+ double stZmax = fpDetInterface->GetZ(iSt) + 0.5 * fpDetInterface->GetThickness(iSt);
+ double stHmin = -fpDetInterface->GetRmax(iSt);
+ double stHmax = +fpDetInterface->GetRmax(iSt);
+
+ auto* pBox = new TBox(stZmin, stHmin, stZmax, stHmax);
+ pBox->SetLineWidth(contWidth);
+ pBox->SetLineStyle(contStyle);
+ pBox->SetLineColor(contColor);
+ pBox->SetFillStyle(contFill);
+ pBox->Draw("SAME");
+ }
+
+ // ** Occupancy in YZ plane
+ MakeCanvas<CbmQaCanvas>("hit_ypos_vs_zpos", "", 600, 400);
+ fvph_hit_ypos_vs_zpos[nSt]->DrawCopy("colz", "");
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ // Station positions in detector IFS
+ double stZmin = fpDetInterface->GetZ(iSt) - 0.5 * fpDetInterface->GetThickness(iSt);
+ double stZmax = fpDetInterface->GetZ(iSt) + 0.5 * fpDetInterface->GetThickness(iSt);
+ double stHmin = -fpDetInterface->GetRmax(iSt);
+ double stHmax = +fpDetInterface->GetRmax(iSt);
+
+ auto* pBox = new TBox(stZmin, stHmin, stZmax, stHmax);
+ pBox->SetLineWidth(contWidth);
+ pBox->SetLineStyle(contStyle);
+ pBox->SetLineColor(contColor);
+ pBox->SetFillStyle(contFill);
+ pBox->Draw("SAME");
+ }
+
+ // ************
+ // ************
+
+ if (IsMCUsed()) {
+
+ // **********************
+ // ** Point occupancies
+
+ // ** Occupancy in XY plane
+ auto* pc_point_ypos_vs_xpos = MakeCanvas<CbmQaCanvas>("point_ypos_vs_xpos", "", 1600, 800);
+ pc_point_ypos_vs_xpos->Divide2D(nSt);
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ pc_point_ypos_vs_xpos->cd(iSt + 1);
+ fvph_point_ypos_vs_xpos[iSt]->DrawCopy("colz", "");
+
+ // Reference size of the station saved to the detector interface
+ auto* pInnerCircle = new TEllipse(0., 0., fpDetInterface->GetRmin(iSt));
+ pInnerCircle->SetLineWidth(contWidth);
+ pInnerCircle->SetLineStyle(contStyle);
+ pInnerCircle->SetLineColor(contColor);
+ pInnerCircle->SetFillStyle(contFill);
+ pInnerCircle->Draw("SAME");
+
+ auto* pOuterCircle = new TEllipse(0., 0., fpDetInterface->GetRmax(iSt));
+ pOuterCircle->SetLineWidth(contWidth);
+ pOuterCircle->SetLineStyle(contStyle);
+ pOuterCircle->SetLineColor(contColor);
+ pOuterCircle->SetFillStyle(contFill);
+ pOuterCircle->Draw("SAME");
+
+ // Square boarders of station, which are used for the field approximation
+ double stXmin = +fpDetInterface->GetXmax(iSt);
+ double stXmax = -fpDetInterface->GetXmax(iSt);
+ double stYmin = -fpDetInterface->GetYmax(iSt);
+ double stYmax = +fpDetInterface->GetYmax(iSt);
+
+ auto* pBox = new TBox(stXmin, stYmin, stXmax, stYmax);
+ pBox->SetLineWidth(contWidth);
+ pBox->SetLineStyle(contStyle);
+ pBox->SetLineColor(contColor);
+ pBox->SetFillStyle(contFill);
+ pBox->Draw("SAME");
+ }
+
+ // ** Occupancy in XZ plane
+ MakeCanvas<CbmQaCanvas>("point_xpos_vs_zpos", "", 600, 400);
+ fvph_point_xpos_vs_zpos[nSt]->SetStats(false);
+ fvph_point_xpos_vs_zpos[nSt]->DrawCopy("colz", "");
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ // Station positions in detector IFS
+ double stZmin = fpDetInterface->GetZ(iSt) - 0.5 * fpDetInterface->GetThickness(iSt);
+ double stZmax = fpDetInterface->GetZ(iSt) + 0.5 * fpDetInterface->GetThickness(iSt);
+ double stHmin = -fpDetInterface->GetRmax(iSt);
+ double stHmax = +fpDetInterface->GetRmax(iSt);
+
+ auto* pBox = new TBox(stZmin, stHmin, stZmax, stHmax);
+ pBox->SetLineWidth(contWidth);
+ pBox->SetLineStyle(contStyle);
+ pBox->SetLineColor(contColor);
+ pBox->SetFillStyle(contFill);
+ pBox->Draw("SAME");
+ }
+
+ // ** Occupancy in YZ plane
+ MakeCanvas<CbmQaCanvas>("point_ypos_vs_zpos", "", 600, 400);
+ fvph_point_ypos_vs_zpos[nSt]->SetStats(false);
+ fvph_point_ypos_vs_zpos[nSt]->DrawCopy("colz", "");
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ // Station positions in detector IFS
+ double stZmin = fpDetInterface->GetZ(iSt) - 0.5 * fpDetInterface->GetThickness(iSt);
+ double stZmax = fpDetInterface->GetZ(iSt) + 0.5 * fpDetInterface->GetThickness(iSt);
+ double stHmin = -fpDetInterface->GetRmax(iSt);
+ double stHmax = +fpDetInterface->GetRmax(iSt);
+
+ auto* pBox = new TBox(stZmin, stHmin, stZmax, stHmax);
+ pBox->SetLineWidth(contWidth);
+ pBox->SetLineStyle(contStyle);
+ pBox->SetLineColor(contColor);
+ pBox->SetFillStyle(contFill);
+ pBox->Draw("SAME");
+ }
+
+ // **************
+ // ** Residuals
+
+ // x-coordinate
+ auto* pc_res_x = MakeCanvas<CbmQaCanvas>("res_x", "Residuals for x coordinate", 1600, 800);
+ pc_res_x->Divide2D(nSt);
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ pc_res_x->cd(iSt + 1);
+ fvph_res_x[iSt]->DrawCopy("", "");
+ }
+
+ // y-coordinate
+ auto* pc_res_y = MakeCanvas<CbmQaCanvas>("res_y", "Residuals for y coordinate", 1600, 800);
+ pc_res_y->Divide2D(nSt);
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ pc_res_y->cd(iSt + 1);
+ fvph_res_y[iSt]->DrawCopy("", "");
+ }
+
+ // time
+ auto* pc_res_t = MakeCanvas<CbmQaCanvas>("res_t", "Residuals for time", 1600, 800);
+ pc_res_t->Divide2D(nSt);
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ pc_res_t->cd(iSt + 1);
+ fvph_res_t[iSt]->DrawCopy("", "");
+ }
+
+
+ // **********
+ // ** Pulls
+
+ // x-coordinate
+ auto* pc_pull_x = MakeCanvas<CbmQaCanvas>("pull_x", "Pulls for x coordinate", 1600, 800);
+ pc_pull_x->Divide2D(nSt);
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ pc_pull_x->cd(iSt + 1);
+ fvph_pull_x[iSt]->DrawCopy("", "");
+ }
+
+ // y-coordinate
+ auto* pc_pull_y = MakeCanvas<CbmQaCanvas>("pull_y", "Pulls for y coordinate", 1600, 800);
+ pc_pull_y->Divide2D(nSt);
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ pc_pull_y->cd(iSt + 1);
+ fvph_pull_y[iSt]->DrawCopy("", "");
+ }
+
+ // time
+ auto* pc_pull_t = MakeCanvas<CbmQaCanvas>("pull_t", "Pulls for time", 1600, 800);
+ pc_pull_t->Divide2D(nSt);
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ pc_pull_t->cd(iSt + 1);
+ fvph_pull_t[iSt]->DrawCopy("", "");
+ }
+
+
+ // ************************************
+ // ** Hit reconstruction efficiencies
+
+ auto* pc_reco_eff_vs_r =
+ MakeCanvas<CbmQaCanvas>("reco_eff_vs_r", "Hit efficiencies wrt distance from center", 1600, 800);
+ pc_reco_eff_vs_r->Divide2D(nSt);
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ pc_reco_eff_vs_r->cd(iSt + 1);
+ fvpe_reco_eff_vs_r[iSt]->Paint("AP");
+ auto* pGr = dynamic_cast<TGraphAsymmErrors*>(fvpe_reco_eff_vs_r[iSt]->GetPaintedGraph());
+ if (!pGr) {
+ LOG(error) << fName << ": unable to get painted graph from efficiency " << fvpe_reco_eff_vs_xy[iSt]->GetName();
+ continue;
+ }
+ pGr->DrawClone("AP");
+ }
+
+ auto* pc_reco_eff_vs_xy = MakeCanvas<CbmQaCanvas>("reco_eff_vs_xy", "Hit efficiencies wrt x and y", 1600, 800);
+ pc_reco_eff_vs_xy->Divide2D(nSt);
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ pc_reco_eff_vs_xy->cd(iSt + 1);
+ fvpe_reco_eff_vs_xy[iSt]->Paint("colz");
+ auto* pH2 = dynamic_cast<TH2F*>(fvpe_reco_eff_vs_xy[iSt]->GetPaintedHistogram());
+ if (!pH2) {
+ LOG(error) << fName << ": unable to get painted histogram from efficiency "
+ << fvpe_reco_eff_vs_xy[iSt]->GetName();
+ continue;
+ }
+ pH2->DrawCopy("colz", "");
+ }
+ }
+
+ return kSUCCESS;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+InitStatus CbmCaInputQaTof::InitHistograms()
+{
+ int nSt = fpDetInterface->GetNtrackingStations();
+
+ // ----- Histograms initialization
+ fvph_hit_ypos_vs_xpos.resize(nSt + 1, nullptr);
+ fvph_hit_ypos_vs_zpos.resize(nSt + 1, nullptr);
+ fvph_hit_xpos_vs_zpos.resize(nSt + 1, nullptr);
+
+ fvph_hit_station_delta_z.resize(nSt);
+
+ fvph_hit_dx.resize(nSt + 1, nullptr);
+ fvph_hit_dy.resize(nSt + 1, nullptr);
+ fvph_hit_dt.resize(nSt + 1, nullptr);
+
+ for (int iSt = 0; iSt <= nSt; ++iSt) {
+ TString nsuff = (iSt == nSt) ? "" : Form("_st%d", iSt); // Histogram name suffix
+ TString tsuff = (iSt == nSt) ? "" : Form(" in TOF station %d", iSt); // Histogram title suffix
+ TString sN = "";
+ TString sT = "";
+
+ // Hit occupancy
+ sN = (TString) "hit_ypos_vs_xpos" + nsuff;
+ sT = (TString) "Hit occupancy in xy-Plane" + tsuff + ";x_{hit} [cm];y_{hit} [cm]";
+ fvph_hit_ypos_vs_xpos[iSt] = MakeHisto<TH2F>(sN, sT, kNbins, kRHitX[0], kRHitX[1], kNbins, kRHitY[0], kRHitY[1]);
+
+ sN = (TString) "hit_xpos_vs_zpos" + nsuff;
+ sT = (TString) "Hit occupancy in xz-Plane" + tsuff + ";z_{hit} [cm];x_{hit} [cm]";
+ fvph_hit_xpos_vs_zpos[iSt] = MakeHisto<TH2F>(sN, sT, kNbinsZ, kRHitZ[0], kRHitZ[1], kNbins, kRHitX[0], kRHitX[1]);
+
+ sN = (TString) "hit_ypos_vs_zpos" + nsuff;
+ sT = (TString) "Hit occupancy in yz-plane" + tsuff + ";z_{hit} [cm];y_{hit} [cm]";
+ fvph_hit_ypos_vs_zpos[iSt] = MakeHisto<TH2F>(sN, sT, kNbinsZ, kRHitZ[0], kRHitZ[1], kNbins, kRHitY[0], kRHitY[1]);
+
+ // Hit errors
+ sN = (TString) "hit_dx" + nsuff;
+ sT = (TString) "Hit position error along x-axis" + tsuff + ";dx_{hit} [cm]";
+ fvph_hit_dx[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, kRHitDx[0], kRHitDx[1]);
+
+ sN = (TString) "hit_dy" + nsuff;
+ sT = (TString) "Hit position error along y-axis" + tsuff + ";dy_{hit} [cm]";
+ fvph_hit_dy[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, kRHitDy[0], kRHitDy[1]);
+
+ sN = (TString) "hit_dt" + nsuff;
+ sT = (TString) "Hit time error" + tsuff + ";dt_{hit} [ns]";
+ fvph_hit_dt[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, kRHitDt[0], kRHitDt[1]);
+
+ if (iSt == nSt) { continue; } // Following histograms are defined only for separate station
+
+ sN = (TString) "hit_station_delta_z" + nsuff;
+ sT = (TString) "Different between hit and station z-positions" + tsuff + ";z_{hit} - z_{st} [cm]";
+ fvph_hit_station_delta_z[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, -5., 5);
+
+ } // loop over station index: end
+
+ // ----- Initialize histograms, which are use MC-information
+ if (IsMCUsed()) {
+ // Resize histogram vectors
+ fvph_n_points_per_hit.resize(nSt + 1, nullptr);
+ fvph_point_ypos_vs_xpos.resize(nSt + 1, nullptr);
+ fvph_point_xpos_vs_zpos.resize(nSt + 1, nullptr);
+ fvph_point_ypos_vs_zpos.resize(nSt + 1, nullptr);
+ fvph_point_hit_delta_z.resize(nSt + 1, nullptr);
+ fvph_res_x.resize(nSt + 1, nullptr);
+ fvph_res_y.resize(nSt + 1, nullptr);
+ fvph_res_t.resize(nSt + 1, nullptr);
+ fvph_pull_x.resize(nSt + 1, nullptr);
+ fvph_pull_y.resize(nSt + 1, nullptr);
+ fvph_pull_t.resize(nSt + 1, nullptr);
+ fvph_res_x_vs_x.resize(nSt + 1, nullptr);
+ fvph_res_y_vs_y.resize(nSt + 1, nullptr);
+ fvph_res_t_vs_t.resize(nSt + 1, nullptr);
+ fvph_pull_x_vs_x.resize(nSt + 1, nullptr);
+ fvph_pull_y_vs_y.resize(nSt + 1, nullptr);
+ fvph_pull_t_vs_t.resize(nSt + 1, nullptr);
+ fvpe_reco_eff_vs_xy.resize(nSt + 1, nullptr);
+ fvpe_reco_eff_vs_r.resize(nSt + 1, nullptr);
+
+ for (int iSt = 0; iSt <= nSt; ++iSt) {
+ TString nsuff = (iSt == nSt) ? "" : Form("_st%d", iSt); // Histogram name suffix
+ TString tsuff = (iSt == nSt) ? "" : Form(" in TOF station %d", iSt); // Histogram title suffix
+ TString sN = "";
+ TString sT = "";
+
+ sN = (TString) "n_points_per_hit" + nsuff;
+ sT = (TString) "Number of points per hit" + tsuff + ";N_{point}/hit";
+ fvph_n_points_per_hit[iSt] = MakeHisto<TH1F>(sN, sT, 10, -0.5, 9.5);
+
+ // Point occupancy
+ sN = (TString) "point_ypos_vs_xpos" + nsuff;
+ sT = (TString) "Point occupancy in XY plane" + tsuff + ";x_{MC} [cm];y_{MC} [cm]";
+ fvph_point_ypos_vs_xpos[iSt] =
+ MakeHisto<TH2F>(sN, sT, kNbins, kRHitX[0], kRHitX[1], kNbins, kRHitY[0], kRHitY[1]);
+
+ sN = (TString) "point_xpos_vs_zpos" + nsuff;
+ sT = (TString) "Point Occupancy in XZ plane" + tsuff + ";z_{MC} [cm];x_{MC} [cm]";
+ fvph_point_xpos_vs_zpos[iSt] =
+ MakeHisto<TH2F>(sN, sT, kNbinsZ, kRHitZ[0], kRHitZ[1], kNbins, kRHitX[0], kRHitX[1]);
+
+ sN = (TString) "point_ypos_vs_zpos" + nsuff;
+ sT = (TString) "Point Occupancy in YZ Plane" + tsuff + ";z_{MC} [cm];y_{MC} [cm]";
+ fvph_point_ypos_vs_zpos[iSt] =
+ MakeHisto<TH2F>(sN, sT, kNbinsZ, kRHitZ[0], kRHitZ[1], kNbins, kRHitY[0], kRHitY[1]);
+
+ // Difference between MC input z and hit z coordinates
+ sN = (TString) "point_hit_delta_z" + nsuff;
+ sT = (TString) "Distance between point and hit along z axis" + tsuff + ";z_{reco} - z_{MC} [cm]";
+ fvph_point_hit_delta_z[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, -0.04, 0.04);
+
+ sN = (TString) "res_x" + nsuff;
+ sT = (TString) "Residuals for x-coordinate" + tsuff + ";x_{reco} - x_{MC} [cm]";
+ fvph_res_x[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, kRResX[0], kRResX[1]);
+
+ sN = (TString) "res_y" + nsuff;
+ sT = (TString) "Residuals for y-coordinate" + tsuff + ";y_{reco} - y_{MC} [cm]";
+ fvph_res_y[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, kRResY[0], kRResY[1]);
+
+ sN = (TString) "res_t" + nsuff;
+ sT = (TString) "Residuals for time" + tsuff + ";t_{reco} - t_{MC} [ns]";
+ fvph_res_t[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, kRResT[0], kRResT[1]);
+
+ sN = (TString) "pull_x" + nsuff;
+ sT = (TString) "Pulls for x-coordinate" + tsuff + ";(x_{reco} - x_{MC}) / #sigma_{x}^{reco}";
+ fvph_pull_x[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, kRPullX[0], kRPullX[1]);
+
+ sN = (TString) "pull_y" + nsuff;
+ sT = (TString) "Pulls for y-coordinate" + tsuff + ";(y_{reco} - y_{MC}) / #sigma_{y}^{reco}";
+ fvph_pull_y[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, kRPullY[0], kRPullY[1]);
+
+ sN = (TString) "pull_t" + nsuff;
+ sT = (TString) "Pulls for time" + tsuff + ";(t_{reco} - t_{MC}) / #sigma_{t}^{reco}";
+ fvph_pull_t[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, kRPullT[0], kRPullT[1]);
+
+ sN = (TString) "res_x_vs_x" + nsuff;
+ sT = (TString) "Residuals for x-coordinate" + tsuff + ";x_{reco} [cm];x_{reco} - x_{MC} [cm]";
+ fvph_res_x_vs_x[iSt] = MakeHisto<TH2F>(sN, sT, kNbins, 0, 0, kNbins, kRResX[0], kRResX[1]);
+
+ sN = (TString) "res_y_vs_y" + nsuff;
+ sT = (TString) "Residuals for y-coordinate" + tsuff + ";y_{reco} [cm];y_{reco} - y_{MC} [cm]";
+ fvph_res_y_vs_y[iSt] = MakeHisto<TH2F>(sN, sT, kNbins, 0, 0, kNbins, kRResY[0], kRResY[1]);
+
+ sN = (TString) "res_t_vs_t" + nsuff;
+ sT = (TString) "Residuals for time" + tsuff + "t_{reco} [ns];t_{reco} - t_{MC} [ns]";
+ fvph_res_t_vs_t[iSt] = MakeHisto<TH2F>(sN, sT, kNbins, 0, 0, kNbins, kRResT[0], kRResT[1]);
+
+ sN = (TString) "pull_x_vs_x" + nsuff;
+ sT = (TString) "Pulls for x-coordinate" + tsuff + "x_{reco} [cm];(x_{reco} - x_{MC}) / #sigma_{x}^{reco}";
+ fvph_pull_x_vs_x[iSt] = MakeHisto<TH2F>(sN, sT, kNbins, 0, 0, kNbins, kRPullX[0], kRPullX[1]);
+
+ sN = (TString) "pull_y_vs_y" + nsuff;
+ sT = (TString) "Pulls for y-coordinate" + tsuff + ";y_{reco} [cm];(y_{reco} - y_{MC}) / #sigma_{y}^{reco}";
+ fvph_pull_y_vs_y[iSt] = MakeHisto<TH2F>(sN, sT, kNbins, 0, 0, kNbins, kRPullY[0], kRPullY[1]);
+
+ sN = (TString) "pull_t_vs_t" + nsuff;
+ sT = (TString) "Pulls for time" + tsuff + ";t_{reco} [ns];(t_{reco} - t_{MC}) / #sigma_{t}^{reco}";
+ fvph_pull_t_vs_t[iSt] = MakeHisto<TH2F>(sN, sT, kNbins, 0, 0, kNbins, kRPullT[0], kRPullT[1]);
+
+
+ sN = (TString) "reco_eff_vs_xy" + nsuff;
+ sT = (TString) "Hit rec. efficiency" + tsuff + ";x_{MC} [cm];y_{MC} [cm];#epsilon";
+ fvpe_reco_eff_vs_xy[iSt] = MakeEfficiency<CbmQaEff>(sN, sT, 100, -50, 50, 100, -50, 50);
+
+ sN = (TString) "reco_eff_vs_r" + nsuff;
+ sT = (TString) "Hit rec. efficiency" + tsuff + ";r_{MC} [cm];#epsilon";
+ fvpe_reco_eff_vs_r[iSt] = MakeEfficiency<CbmQaEff>(sN, sT, 100, 0, 100);
+ }
+ }
+ return kSUCCESS;
+}
diff --git a/reco/L1/qa/CbmCaInputQaTof.h b/reco/L1/qa/CbmCaInputQaTof.h
new file mode 100644
index 0000000000..3bf035ced9
--- /dev/null
+++ b/reco/L1/qa/CbmCaInputQaTof.h
@@ -0,0 +1,178 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// \file CbmCaInputQaTof.h
+/// \date 30.01.2023
+/// \brief QA-task for CA tracking input from TOF detector (header)
+/// \author S.Zharko <s.zharko@gsi.de>
+
+
+#ifndef CbmCaInputQaTof_h
+#define CbmCaInputQaTof_h 1
+
+#include "CbmCaInputQaBase.h"
+#include "CbmMCDataManager.h"
+#include "CbmQaTask.h"
+
+#include "TMath.h"
+
+#include <set>
+#include <unordered_map>
+#include <vector>
+
+class CbmMCEventList;
+class CbmMCDataArray;
+class CbmMCDataManager;
+class CbmTimeSlice;
+class CbmMatch;
+class CbmTofHit;
+class CbmTofPoint;
+class CbmTofTrackingInterface;
+class TClonesArray;
+class TH1F;
+class CbmQaEff;
+
+
+/// A QA-task class, which provides assurance of TOF hits and geometry
+///
+class CbmCaInputQaTof : public CbmQaTask, public CbmCaInputQaBase {
+public:
+ /// Constructor from parameters
+ /// \param verbose Verbose level
+ /// \param isMCUsed Flag, whether MC information is available for this task
+ CbmCaInputQaTof(int verbose, bool isMCUsed);
+
+ ClassDef(CbmCaInputQaTof, 0);
+
+protected:
+ // ********************************************
+ // ** Virtual method override from CbmQaTask **
+ // ********************************************
+
+ /// Checks results of the QA
+ /// \return Success flag
+ bool Check();
+
+ /// Initializes data branches
+ InitStatus InitDataBranches();
+
+ /// Initializes canvases
+ InitStatus InitCanvases();
+
+ /// Initializes histograms
+ InitStatus InitHistograms();
+
+ /// Fills histograms per event or time-slice
+ void FillHistograms();
+
+ /// De-initializes histograms
+ void DeInit();
+
+private:
+ // *********************
+ // ** Private methods **
+ // *********************
+
+ /// Gets index of station by pointer to MC point
+ /// \param pPoint Pointer to TOF MC point
+ int GetStationID(const CbmTofPoint* pPoint) const;
+
+ /// Fill vector of interactions and corresponding match objects
+ /// \param pvInters Pointer to interactions array
+ /// \param pvMatches Pointer to hit->interaction matches
+ void FillInteractions(std::shared_ptr<TClonesArray>& pvInters, std::shared_ptr<TClonesArray>& pvMatches);
+
+ // *********************
+ // ** Class variables **
+ // *********************
+
+
+ // ** Data branches **
+ CbmTofTrackingInterface* fpDetInterface = nullptr; ///< Instance of detector interface
+
+ CbmTimeSlice* fpTimeSlice = nullptr; ///< Pointer to current time-slice
+
+ TClonesArray* fpHits = nullptr; ///< Array of hits
+
+ CbmMCDataManager* fpMCDataManager = nullptr; ///< MC data manager
+ CbmMCEventList* fpMCEventList = nullptr; ///< MC event list
+ CbmMCDataArray* fpMCTracks = nullptr; ///< Array of MC tracks
+ CbmMCDataArray* fpMCPoints = nullptr; ///< Array of MC points
+
+ TClonesArray* fpHitMatches = nullptr; ///< Array of hit matches
+
+ // ** Histograms **
+
+ static constexpr double kEffRangeMin = 10.; ///< Lower limit of hit distance for efficiency integration [cm]
+ static constexpr double kEffRangeMax = 30.; ///< Upper limit of hit distance for efficiency integration [cm]
+
+ static constexpr double kRHitX[2] = {-200., 200.}; ///< Range for hit x coordinate [cm]
+ static constexpr double kRHitY[2] = {-200., 200.}; ///< Range for hit y coordinate [cm]
+ static constexpr double kRHitZ[2] = {750., 850.}; ///< Range for hit z coordinate [cm]
+
+ static constexpr int kNbins = 200; ///< General number of bins
+ static constexpr int kNbinsZ = 800; ///< Number of bins for z coordinate axis
+
+ static constexpr double kRHitDx[2] = {-.05, .005}; ///< Range for hit x coordinate [cm]
+ static constexpr double kRHitDy[2] = {-.05, .005}; ///< Range for hit y coordinate [cm]
+ static constexpr double kRHitDt[2] = {-10., 10.}; ///< Range for hit time [ns]
+
+ static constexpr double kRResX[2] = {-2., 2.}; ///< Range for residual of x coordinate [cm]
+ static constexpr double kRResY[2] = {-2., 2.}; ///< Range for residual of y coordinate [cm]
+ static constexpr double kRResT[2] = {-.5, .5}; ///< Range for residual of time [ns]
+
+ static constexpr double kRPullX[2] = {-10., 10.}; ///< Range for pull of x coordinate
+ static constexpr double kRPullY[2] = {-10., 10.}; ///< Range for pull of y coordinate
+ static constexpr double kRPullT[2] = {-5., 5.}; ///< Range for pull of time
+
+
+ // NOTE: the last element of each vector stands for integral distribution over all stations
+
+ // hit occupancy
+ std::vector<TH2F*> fvph_hit_ypos_vs_xpos; ///< hit ypos vs xpos in different stations
+ std::vector<TH2F*> fvph_hit_xpos_vs_zpos; ///< hit xpos vs zpos in different stations
+ std::vector<TH2F*> fvph_hit_ypos_vs_zpos; ///< hit ypos vs zpos in different stations
+
+ // difference between hit and station z positions
+ std::vector<TH1F*> fvph_hit_station_delta_z; ///< Difference between station and hit z positions [cm]
+
+ // hit errors
+ std::vector<TH1F*> fvph_hit_dx;
+ std::vector<TH1F*> fvph_hit_dy;
+ std::vector<TH1F*> fvph_hit_dt;
+
+ // MC points occupancy
+ std::vector<TH1F*> fvph_n_points_per_hit; ///< number of points per one hit in different stations
+
+
+ std::vector<TH2F*> fvph_point_ypos_vs_xpos; ///< point ypos vs xpos in different stations
+ std::vector<TH2F*> fvph_point_xpos_vs_zpos; ///< point xpos vs zpos in different stations
+ std::vector<TH2F*> fvph_point_ypos_vs_zpos; ///< point ypos vs zpos in different stations
+
+ std::vector<TH1F*> fvph_point_hit_delta_z; ///< difference between z of hit and MC point in different stations
+
+ // Residuals
+ std::vector<TH1F*> fvph_res_x; ///< residual for x coordinate in different stations
+ std::vector<TH1F*> fvph_res_y; ///< residual for y coordinate in different stations
+ std::vector<TH1F*> fvph_res_t; ///< residual for t coordinate in different stations
+
+ std::vector<TH2F*> fvph_res_x_vs_x; ///< residual for x coordinate in different stations
+ std::vector<TH2F*> fvph_res_y_vs_y; ///< residual for y coordinate in different stations
+ std::vector<TH2F*> fvph_res_t_vs_t; ///< residual for t coordinate in different stations
+
+ // Pulls
+ std::vector<TH1F*> fvph_pull_x; ///< pull for x coordinate in different stations
+ std::vector<TH1F*> fvph_pull_y; ///< pull for y coordinate in different stations
+ std::vector<TH1F*> fvph_pull_t; ///< pull for t coordinate in different stations
+
+ std::vector<TH2F*> fvph_pull_x_vs_x; ///< pull for x coordinate in different stations
+ std::vector<TH2F*> fvph_pull_y_vs_y; ///< pull for y coordinate in different stations
+ std::vector<TH2F*> fvph_pull_t_vs_t; ///< pull for t coordinate in different stations
+
+ // Hit efficiencies
+ std::vector<CbmQaEff*> fvpe_reco_eff_vs_xy; ///< Efficiency of hit reconstruction vs. x and y coordinates of MC
+ std::vector<CbmQaEff*> fvpe_reco_eff_vs_r; ///< Efficiency of hit reconstruction vs. distance from center
+};
+
+#endif // CbmCaInputQaTof_h
diff --git a/reco/L1/qa/CbmCaInputQaTrd.cxx b/reco/L1/qa/CbmCaInputQaTrd.cxx
new file mode 100644
index 0000000000..e97fa1a4a1
--- /dev/null
+++ b/reco/L1/qa/CbmCaInputQaTrd.cxx
@@ -0,0 +1,1056 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// \file CbmCaInputQaTrd.cxx
+/// \date 13.01.2023
+/// \brief QA-task for CA tracking input from TRD detector (implementation)
+/// \author S.Zharko <s.zharko@gsi.de>
+
+#include "CbmCaInputQaTrd.h"
+
+#include "CbmAddress.h"
+#include "CbmMCDataArray.h"
+#include "CbmMCEventList.h"
+#include "CbmMCTrack.h"
+#include "CbmMatch.h"
+#include "CbmQaCanvas.h"
+#include "CbmQaEff.h"
+#include "CbmTimeSlice.h"
+#include "CbmTrdHit.h"
+#include "CbmTrdPoint.h"
+#include "CbmTrdTrackingInterface.h"
+
+#include "FairRootManager.h"
+#include "Logger.h"
+
+#include "TBox.h"
+#include "TClonesArray.h"
+#include "TEllipse.h"
+#include "TF1.h"
+#include "TFormula.h"
+#include "TGraphAsymmErrors.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TMath.h"
+#include "TStyle.h"
+
+#include <algorithm>
+#include <fstream>
+#include <iomanip>
+#include <numeric>
+
+#include "L1Constants.h"
+
+ClassImp(CbmCaInputQaTrd);
+
+namespace phys = L1Constants::phys; // from L1Constants.h
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+CbmCaInputQaTrd::CbmCaInputQaTrd(int verbose, bool isMCUsed) : CbmQaTask("CbmCaInputQaTrd", "catrd", verbose, isMCUsed)
+{
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+bool CbmCaInputQaTrd::Check()
+{
+ bool res = true;
+
+ int nSt = fpDetInterface->GetNtrackingStations();
+
+
+ // **************************************************************
+ // ** Basic checks, available both for real and simulated data **
+ // **************************************************************
+
+ // ----- Checks for mismatches in the ordering of the stations
+ //
+ std::vector<double> vStationPos(nSt, 0.);
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ vStationPos[iSt] = fpDetInterface->GetZ(iSt);
+ }
+
+ if (!std::is_sorted(vStationPos.cbegin(), vStationPos.cend(), [](int l, int r) { return l <= r; })) {
+ if (fVerbose > 0) {
+ LOG(error) << fName << ": stations are ordered improperly along the beam axis:";
+ for (auto zPos : vStationPos) {
+ LOG(error) << "\t- " << zPos;
+ }
+ }
+ res = false;
+ }
+
+ // ----- Checks for mismatch between station and hit z positions
+ // The purpose of this block is to be ensured, that hits belong to the correct tracking station. For each tracking
+ // station a unified position z-coordinate is defined, which generally differs from the corresponding positions of
+ // reconstructed hits. This is due to non-regular position along the beam axis for each detector sensor. Thus, the
+ // positions of the hits along z-axis are distributed relatively to the position of the station in some range.
+ // If hits goes out this range, it can signal about a mismatch between hits and geometry. For limits of the range
+ // one should select large enough values to cover the difference described above and in the same time small enough
+ // to avoid overlaps with the neighboring stations.
+ // For each station, a distribution of z_{hit} - z_{st} is integrated over the defined range and scaled by the
+ // total number of entries to the distribution. If this value is smaller then unity, then some of the hits belong to
+ // another station.
+ //
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ int nHits = fvph_hit_station_delta_z[iSt]->GetEntries();
+ if (!nHits) {
+ LOG_IF(error, fVerbose > 0) << fName << ": station " << iSt << " does not have hits";
+ res = false;
+ continue;
+ }
+ int iBinMin = fvph_hit_station_delta_z[iSt]->FindBin(-fMaxDiffZStHit);
+ int iBinMax = fvph_hit_station_delta_z[iSt]->FindBin(+fMaxDiffZStHit);
+
+ if (fvph_hit_station_delta_z[iSt]->Integral(iBinMin, iBinMax) < nHits) {
+ LOG_IF(error, fVerbose > 0) << fName << ": station " << iSt << " has mismatches in hit z-positions";
+ res = false;
+ }
+ }
+
+
+ // *******************************************************
+ // ** Additional checks, if MC information is available **
+ // *******************************************************
+
+ if (IsMCUsed()) {
+ // ----- Check efficiencies
+ // Error is raised, if any station has integrated efficiency lower then a selected threshold.
+ // Integrated efficiency is estimated with fitting the efficiency curve for hit distances (r) with a uniform
+ // distribution in the range from kEffRangeMin to kEffRangeMax, where the efficiency is assigned to be constant
+ //
+ // Fit efficiency curves
+ LOG(info) << "-- Hit efficiency integrated over hit distance from station center";
+ LOG(info) << "\tintegration range: [" << fEffRange[0] << ", " << fEffRange[1] << "] cm";
+
+ auto* pEffTable = MakeTable("eff_table", "Efficiency table", nSt, 2);
+ pEffTable->SetNamesOfCols({"Station ID", "Efficiency"});
+ pEffTable->SetColWidth(20);
+
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ auto pEff = std::unique_ptr<CbmQaEff>(fvpe_reco_eff_vs_r[iSt]->Integral(fEffRange[0], fEffRange[1]));
+ double eff = pEff->GetEfficiency(1);
+ pEffTable->SetCell(iSt, 0, iSt);
+ pEffTable->SetCell(iSt, 1, eff);
+ res = CheckRange("Hit finder efficiency in station " + std::to_string(iSt), eff, fEffThrsh, 1.000);
+ }
+
+ LOG(info) << '\n' << pEffTable->ToString(3);
+
+ // ----- Checks for residuals
+ // Check hits for potential biases from central values
+
+ // Function to fit a residual distribution, returns a structure
+ auto FitResidualsAndGetMean = [&](TH1* pHist) {
+ auto* pFit = (TF1*) gROOT->FindObject("gausn");
+ LOG_IF(fatal, !pFit) << fName << ": residuals fit function is not found";
+ double statMean = pHist->GetMean();
+ double statSigm = pHist->GetStdDev();
+ pFit->SetParameters(100., statMean, statSigm);
+ pHist->Fit(pFit, "MQ");
+ pHist->Fit(pFit, "MQ");
+ pHist->Fit(pFit, "MQ");
+ // NOTE: Several fit procedures are made to avoid empty fit results
+ std::array<double, 3> result;
+ result[0] = pFit->GetParameter(1);
+ result[1] = -pFit->GetParameter(2) * fResMeanThrsh;
+ result[2] = +pFit->GetParameter(2) * fResMeanThrsh;
+ return result;
+ };
+
+ auto* pResidualsTable = MakeTable("residuals_mean", "Residual mean values in different stations", nSt, 4);
+ pResidualsTable->SetNamesOfCols({"Station ID", "Residual(x) [cm]", "Residual(y) [cm]", "Residual(t) [ns]"});
+ pResidualsTable->SetColWidth(20);
+
+ // Fill residuals fit results and check boundaries
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ auto fitX = FitResidualsAndGetMean(fvph_res_x[iSt]);
+ auto fitY = FitResidualsAndGetMean(fvph_res_y[iSt]);
+ auto fitT = FitResidualsAndGetMean(fvph_res_t[iSt]);
+ res = CheckRange("Mean of x residual [cm] in station " + std::to_string(iSt), fitX[0], fitX[1], fitX[2]) && res;
+ res = CheckRange("Mean of y residual [cm] in station " + std::to_string(iSt), fitY[0], fitY[1], fitY[2]) && res;
+ res = CheckRange("Mean of t residual [ns] in station " + std::to_string(iSt), fitT[0], fitT[1], fitT[2]) && res;
+ pResidualsTable->SetCell(iSt, 0, iSt);
+ pResidualsTable->SetCell(iSt, 1, fitX[0]);
+ pResidualsTable->SetCell(iSt, 2, fitX[1]);
+ pResidualsTable->SetCell(iSt, 3, fitX[2]);
+ }
+
+ LOG(info) << '\n' << pResidualsTable->ToString(8);
+
+ // ----- Checks for pulls
+ // For the hit pull is defined as a ration of the difference between hit and MC-point position or time component
+ // values and an error of the hit position (or time) component, calculated in the same z-positions. In ideal case,
+ // when the resolutions are well determined for detecting stations, the pull distributions should have a RMS equal
+ // to unity. Here we allow a RMS of the pull distributions to be varied in a predefined range. If the RMS runs out
+ // this range, QA task fails.
+ // TODO: Add checks for center
+
+ // Fit pull distributions
+
+ // ** Kaniadakis Gaussian distribution, gives smaller chi2 / NDF
+ //if (!gROOT->FindObject("Expk")) {
+ // new TFormula("Expk", "TMath::Power(TMath::Sqrt(1 + x[1] * x[1] * x[0] * x[0]) + x[0] * x[1], 1./x[1])");
+ //}
+ //TF1* pPullFit = new TF1("pullFitGausn", "[0] * Expk(-[2] * (x - [1]) * (x - [1]), [3])", -10., 10.);
+
+ auto FitPullsAndGetSigma = [&](TH1* pHist) {
+ auto* pFit = (TF1*) gROOT->FindObject("gausn");
+ LOG_IF(fatal, !pFit) << fName << ": pulls fit function is not found";
+ pFit->SetParameters(100, 0.001, 1.000);
+ pHist->Fit(pFit, "Q");
+ return pFit->GetParameter(2);
+ };
+
+ auto* pPullsTable = MakeTable("pulls_rms", "Pulls std. dev. values in different stations", nSt, 4);
+ pPullsTable->SetNamesOfCols({"Station ID", "Pull(x) sigm", "Pull(y) sigm", "Pull(z) sigm"});
+ pPullsTable->SetColWidth(20);
+
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ double rmsPullX = FitPullsAndGetSigma(fvph_pull_x[iSt]);
+ double rmsPullY = FitPullsAndGetSigma(fvph_pull_y[iSt]);
+ double rmsPullT = FitPullsAndGetSigma(fvph_pull_t[iSt]);
+
+ double rmsPullHi = 1. + fPullWidthThrsh;
+ double rmsPullLo = 1. - fPullWidthThrsh;
+
+ res = CheckRange("Rms for x pull in station " + std::to_string(iSt), rmsPullX, rmsPullLo, rmsPullHi) && res;
+ res = CheckRange("Rms for y pull in station " + std::to_string(iSt), rmsPullY, rmsPullLo, rmsPullHi) && res;
+ res = CheckRange("Rms for t pull in station " + std::to_string(iSt), rmsPullT, rmsPullLo, rmsPullHi) && res;
+ pPullsTable->SetCell(iSt, 0, iSt);
+ pPullsTable->SetCell(iSt, 1, rmsPullX);
+ pPullsTable->SetCell(iSt, 2, rmsPullY);
+ pPullsTable->SetCell(iSt, 3, rmsPullT);
+ }
+
+ LOG(info) << '\n' << pPullsTable->ToString(3);
+ }
+
+ return res;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void CbmCaInputQaTrd::DeInit()
+{
+ // Vectors with pointers to histograms
+ fvph_hit_ypos_vs_xpos.clear();
+ fvph_hit_xpos_vs_zpos.clear();
+ fvph_hit_ypos_vs_zpos.clear();
+
+ fvph_hit_station_delta_z.clear();
+
+ fvph_hit_dx.clear();
+ fvph_hit_dy.clear();
+ fvph_hit_dt.clear();
+
+ fvph_n_points_per_hit.clear();
+
+ fvph_point_ypos_vs_xpos.clear();
+ fvph_point_xpos_vs_zpos.clear();
+ fvph_point_ypos_vs_zpos.clear();
+
+ fvph_point_hit_delta_z.clear();
+
+ fvph_res_x.clear();
+ fvph_res_y.clear();
+ fvph_res_t.clear();
+
+ fvph_pull_x.clear();
+ fvph_pull_y.clear();
+ fvph_pull_t.clear();
+
+ fvph_res_x_vs_x.clear();
+ fvph_res_y_vs_y.clear();
+ fvph_res_t_vs_t.clear();
+
+ fvph_pull_x_vs_x.clear();
+ fvph_pull_y_vs_y.clear();
+ fvph_pull_t_vs_t.clear();
+
+ fvpe_reco_eff_vs_xy.clear();
+ fvpe_reco_eff_vs_r.clear();
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void CbmCaInputQaTrd::FillHistograms()
+{
+ int nSt = fpDetInterface->GetNtrackingStations();
+ int nHits = fpHits->GetEntriesFast();
+ int nMCevents = (IsMCUsed()) ? fpMCEventList->GetNofEvents() : -1;
+
+ std::vector<std::vector<int>> vNofHitsPerPoint; // Number of hits per MC point in different MC events
+
+ if (IsMCUsed()) {
+ vNofHitsPerPoint.resize(nMCevents);
+ for (int iE = 0; iE < nMCevents; ++iE) {
+ int iFile = fpMCEventList->GetFileIdByIndex(iE);
+ int iEvent = fpMCEventList->GetEventIdByIndex(iE);
+ int nPoints = fpMCPoints->Size(iFile, iEvent);
+ vNofHitsPerPoint[iE].resize(nPoints, 0);
+ }
+ }
+
+ for (int iH = 0; iH < nHits; ++iH) {
+ const auto* pHit = dynamic_cast<const CbmTrdHit*>(fpHits->At(iH));
+ LOG_IF(fatal, !pHit) << fName << ": hit with iH = " << iH << " is not an CbmTrdHit (dynamic cast failed)";
+
+ // *************************
+ // ** Reconstructed hit QA
+
+ // Hit station geometry info
+ int iSt = fpDetInterface->GetTrackingStationIndex(pHit);
+ LOG_IF(fatal, iSt < 0 || iSt >= nSt) << fName << ": index of station (" << iSt << ") is out of range "
+ << "for hit with id = " << iH;
+
+ int hitType = pHit->GetClassType(); // TRD-1D, TRD-2D
+
+ // Hit position
+ double xHit = pHit->GetX();
+ double yHit = pHit->GetY();
+ double zHit = pHit->GetZ();
+
+ double dxHit = pHit->GetDx();
+ double dyHit = pHit->GetDy();
+ //double dxyHit = pHit->GetDxy();
+
+ // Hit time
+ double tHit = pHit->GetTime();
+ double dtHit = pHit->GetTimeError();
+
+ fvph_hit_ypos_vs_xpos[iSt]->Fill(xHit, yHit);
+ fvph_hit_xpos_vs_zpos[iSt]->Fill(zHit, xHit);
+ fvph_hit_ypos_vs_zpos[iSt]->Fill(zHit, yHit);
+
+ fvph_hit_station_delta_z[iSt]->Fill(zHit - fpDetInterface->GetZ(iSt));
+
+ fvph_hit_dx[iSt]->Fill(dxHit);
+ fvph_hit_dy[iSt]->Fill(dyHit);
+ fvph_hit_dt[iSt]->Fill(dtHit);
+
+ fvph_hit_ypos_vs_xpos[nSt + hitType]->Fill(xHit, yHit);
+ fvph_hit_xpos_vs_zpos[nSt + hitType]->Fill(zHit, xHit);
+ fvph_hit_ypos_vs_zpos[nSt + hitType]->Fill(zHit, yHit);
+ fvph_hit_dx[nSt + hitType]->Fill(dxHit);
+ fvph_hit_dy[nSt + hitType]->Fill(dyHit);
+ fvph_hit_dt[nSt + hitType]->Fill(dtHit);
+
+
+ // **********************
+ // ** MC information QA
+
+ if (IsMCUsed()) {
+ CbmMatch* pHitMatch = dynamic_cast<CbmMatch*>(fpHitMatches->At(iH));
+ LOG_IF(fatal, !pHitMatch) << fName << ": match object not found for hit ID " << iH;
+
+ // Evaluate number of hits per one MC point
+ int nMCpoints = 0; // Number of MC points for this hit
+
+ for (int iLink = 0; iLink < pHitMatch->GetNofLinks(); ++iLink) {
+ const auto& link = pHitMatch->GetLink(iLink);
+
+ int iP = link.GetIndex(); // Index of MC point
+
+ // Skip noisy links
+ if (iP < 0) { continue; }
+
+ ++nMCpoints;
+
+ int iE = fpMCEventList->GetEventIndex(link);
+
+ LOG_IF(fatal, iE < 0 || iE >= nMCevents) << fName << ": id of MC event is out of range (hit id = " << iH
+ << ", link id = " << iLink << ", event id = " << iE << ')';
+
+ LOG_IF(fatal, iP >= (int) vNofHitsPerPoint[iE].size())
+ << fName << ": index of MC point is out of range (hit id = " << iH << ", link id = " << iLink
+ << ", event id = " << iE << ", point id = " << iP << ')';
+ vNofHitsPerPoint[iE][iP]++;
+ }
+
+ fvph_n_points_per_hit[iSt]->Fill(nMCpoints);
+ fvph_n_points_per_hit[nSt + hitType]->Fill(nMCpoints);
+
+ if (nMCpoints != 1) { continue; } // ??
+
+ // The best link to in the match (probably, the cut on nMCpoints is meaningless)
+ const auto& bestPointLink = pHitMatch->GetMatchedLink();
+
+ // Skip noisy links
+ if (bestPointLink.GetIndex() < 0) { continue; }
+
+ // Point matched by the best link
+ const auto* pMCPoint = dynamic_cast<const CbmTrdPoint*>(fpMCPoints->Get(bestPointLink));
+ LOG_IF(fatal, !pMCPoint) << fName << ": MC point object does not exist for hit " << iH;
+
+ // MC track for this point
+ CbmLink bestTrackLink = bestPointLink;
+ bestTrackLink.SetIndex(pMCPoint->GetTrackID());
+ const auto* pMCTrack = dynamic_cast<const CbmMCTrack*>(fpMCTracks->Get(bestTrackLink));
+ LOG_IF(fatal, !pMCTrack) << fName << ": MC track object does not exist for hit " << iH << " and link: ";
+
+ double t0MC = fpMCEventList->GetEventTime(bestPointLink);
+ LOG_IF(fatal, t0MC < 0) << fName << ": MC time zero is lower then 0 ns: " << t0MC;
+
+
+ // ----- MC point properties
+ //
+ double mass = pMCTrack->GetMass();
+ //double charge = pMCTrack->GetCharge();
+ //double pdg = pMCTrack->GetPdgCode();
+
+ // Entrance position and time
+ double xMC = pMCPoint->GetX();
+ double yMC = pMCPoint->GetY();
+ double zMC = pMCPoint->GetZ();
+ double tMC = pMCPoint->GetTime() + t0MC;
+
+ // MC point entrance momenta
+ double pxMC = pMCPoint->GetPx();
+ double pyMC = pMCPoint->GetPy();
+ double pzMC = pMCPoint->GetPz();
+ double pMC = sqrt(pxMC * pxMC + pyMC * pyMC + pzMC * pzMC);
+
+ // MC point exit momenta
+ // double pxMCo = pMCPoint->GetPxOut();
+ // double pyMCo = pMCPoint->GetPyOut();
+ // double pzMCo = pMCPoint->GetPzOut();
+ // double pMCo = sqrt(pxMCo * pxMCo + pyMCo * pyMCo + pzMCo * pzMCo);
+
+ // Position and time shifted to z-coordinate of the hit
+ double shiftZ = zHit - zMC; // Difference between hit and point z positions
+ double xMCs = xMC + shiftZ * pxMC / pzMC;
+ double yMCs = yMC + shiftZ * pyMC / pzMC;
+ double tMCs = tMC + shiftZ / (pzMC * phys::kSpeedOfLight) * sqrt(mass * mass + pMC * pMC);
+
+ // Residuals
+ double xRes = xHit - xMCs;
+ double yRes = yHit - yMCs;
+ double tRes = tHit - tMCs;
+
+ // ------ Cuts
+
+ if (std::fabs(pMCPoint->GetPzOut()) < fMinMomentum) { continue; } // CUT ON MINIMUM MOMENTUM
+ //if (pMCo < cuts::kMinP) { continue; } // Momentum threshold
+
+
+ fvph_point_ypos_vs_xpos[iSt]->Fill(xMC, yMC);
+ fvph_point_xpos_vs_zpos[iSt]->Fill(zMC, xMC);
+ fvph_point_ypos_vs_zpos[iSt]->Fill(zMC, yMC);
+
+ fvph_point_hit_delta_z[iSt]->Fill(shiftZ);
+
+ fvph_res_x[iSt]->Fill(xRes);
+ fvph_res_y[iSt]->Fill(yRes);
+ fvph_res_t[iSt]->Fill(tRes);
+
+ fvph_pull_x[iSt]->Fill(xRes / dxHit);
+ fvph_pull_y[iSt]->Fill(yRes / dyHit);
+ fvph_pull_t[iSt]->Fill(tRes / dtHit);
+
+ fvph_res_x_vs_x[iSt]->Fill(xHit, xRes);
+ fvph_res_y_vs_y[iSt]->Fill(yHit, yRes);
+ fvph_res_t_vs_t[iSt]->Fill(tHit, tRes);
+
+ fvph_pull_x_vs_x[iSt]->Fill(xHit, xRes / dxHit);
+ fvph_pull_y_vs_y[iSt]->Fill(yHit, yRes / dyHit);
+ fvph_pull_t_vs_t[iSt]->Fill(tHit, tRes / dtHit);
+
+ fvph_point_ypos_vs_xpos[nSt + hitType]->Fill(xMC, yMC);
+ fvph_point_xpos_vs_zpos[nSt + hitType]->Fill(zMC, xMC);
+ fvph_point_ypos_vs_zpos[nSt + hitType]->Fill(zMC, yMC);
+
+ fvph_point_hit_delta_z[nSt + hitType]->Fill(shiftZ);
+
+ fvph_res_x[nSt + hitType]->Fill(xRes);
+ fvph_res_y[nSt + hitType]->Fill(yRes);
+ fvph_res_t[nSt + hitType]->Fill(tRes);
+
+ fvph_pull_x[nSt + hitType]->Fill(xRes / dxHit);
+ fvph_pull_y[nSt + hitType]->Fill(yRes / dyHit);
+ fvph_pull_t[nSt + hitType]->Fill(tRes / dtHit);
+
+ fvph_res_x_vs_x[nSt + hitType]->Fill(xHit, xRes);
+ fvph_res_y_vs_y[nSt + hitType]->Fill(yHit, yRes);
+ fvph_res_t_vs_t[nSt + hitType]->Fill(tHit, tRes);
+
+ fvph_pull_x_vs_x[nSt + hitType]->Fill(xHit, xRes / dxHit);
+ fvph_pull_y_vs_y[nSt + hitType]->Fill(yHit, yRes / dyHit);
+ fvph_pull_t_vs_t[nSt + hitType]->Fill(tHit, tRes / dtHit);
+ }
+ } // loop over hits: end
+
+ // Fill hit reconstruction efficiencies
+ if (IsMCUsed()) {
+ for (int iE = 0; iE < nMCevents; ++iE) {
+ int iFile = fpMCEventList->GetFileIdByIndex(iE);
+ int iEvent = fpMCEventList->GetEventIdByIndex(iE);
+ int nPoints = fpMCPoints->Size(iFile, iEvent);
+
+ for (int iP = 0; iP < nPoints; ++iP) {
+ const auto* pMCPoint = dynamic_cast<const CbmTrdPoint*>(fpMCPoints->Get(iFile, iEvent, iP));
+ LOG_IF(fatal, !pMCPoint) << fName << ": MC point does not exist for iFile = " << iFile
+ << ", iEvent = " << iEvent << ", iP = " << iP;
+ int address = pMCPoint->GetDetectorID();
+ int iSt = fpDetInterface->GetTrackingStationIndex(address);
+ LOG_IF(fatal, iSt < 0 || iSt >= nSt)
+ << fName << ": MC point for FEI = " << iFile << ", " << iEvent << ", " << iP << " and address " << address
+ << " has wrong station index: iSt = " << iSt;
+
+ double xMC = pMCPoint->GetXIn();
+ double yMC = pMCPoint->GetYIn();
+ double rMC = sqrt(xMC * xMC + yMC * yMC);
+
+ // Conditions under which point is accounted as reconstructed: point
+ bool ifPointHasHits = (vNofHitsPerPoint[iE][iP] > 0);
+
+ fvpe_reco_eff_vs_xy[iSt]->Fill(ifPointHasHits, xMC, yMC);
+ fvpe_reco_eff_vs_xy[nSt]->Fill(ifPointHasHits, xMC, yMC);
+
+ fvpe_reco_eff_vs_r[iSt]->Fill(ifPointHasHits, rMC);
+ fvpe_reco_eff_vs_r[nSt]->Fill(ifPointHasHits, rMC);
+
+ } // loop over MC-points: end
+ } // loop over MC-events: end
+ } // MC is used: end
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+InitStatus CbmCaInputQaTrd::InitDataBranches()
+{
+ // STS tracking detector interface
+ fpDetInterface = CbmTrdTrackingInterface::Instance();
+
+ LOG_IF(fatal, !fpDetInterface) << "\033[1;31m" << fName << ": tracking detector interface is undefined\033[0m";
+
+ // FairRootManager
+ auto* pFairRootManager = FairRootManager::Instance();
+ LOG_IF(fatal, !pFairRootManager) << "\033[1;31m" << fName << ": FairRootManager instance is a null pointer\033[0m";
+
+ // Time-slice
+ fpTimeSlice = dynamic_cast<CbmTimeSlice*>(pFairRootManager->GetObject("TimeSlice."));
+ LOG_IF(fatal, !fpTimeSlice) << "\033[1;31m" << fName << ": time-slice branch is not found\033[0m";
+
+ // ----- Reconstructed data-branches initialization
+
+ // Hits container
+ fpHits = dynamic_cast<TClonesArray*>(pFairRootManager->GetObject("TrdHit"));
+ LOG_IF(fatal, !fpHits) << "\033[1;31m" << fName << ": container of reconstructed hits in TRD is not found\033[0m";
+
+ // ----- MC-information branches initialization
+ if (IsMCUsed()) {
+ // MC manager
+ fpMCDataManager = dynamic_cast<CbmMCDataManager*>(pFairRootManager->GetObject("MCDataManager"));
+ LOG_IF(fatal, !fpMCDataManager) << "\033[1;31m" << fName << ": MC data manager branch is not found\033[0m";
+
+ // MC event list
+ fpMCEventList = dynamic_cast<CbmMCEventList*>(pFairRootManager->GetObject("MCEventList."));
+ LOG_IF(fatal, !fpMCEventList) << "\033[1;31m" << fName << ": MC event list branch is not found\033[0m";
+
+ // MC tracks
+ fpMCTracks = fpMCDataManager->InitBranch("MCTrack");
+ LOG_IF(fatal, !fpMCTracks) << "\033[1;31m" << fName << ": MC track branch is not found\033[0m";
+
+ // MC points
+ fpMCPoints = fpMCDataManager->InitBranch("TrdPoint");
+ LOG_IF(fatal, !fpMCTracks) << "\033[1;31m" << fName << ": MC point branch is not found for TRD\033[0m";
+
+ // Hit matches
+ fpHitMatches = dynamic_cast<TClonesArray*>(pFairRootManager->GetObject("TrdHitMatch"));
+ LOG_IF(fatal, !fpHitMatches) << "\033[1;31m]" << fName << ": hit match branch is not found for TRD\033[0m";
+ }
+
+ return kSUCCESS;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+InitStatus CbmCaInputQaTrd::InitCanvases()
+{
+ gStyle->SetOptFit(1);
+ int nSt = fpDetInterface->GetNtrackingStations();
+
+ // ********************
+ // ** Drawing options
+
+ // Contours
+ constexpr auto contColor = kOrange + 7;
+ constexpr auto contWidth = 2; // Line width
+ constexpr auto contStyle = 2; // Line style
+ constexpr auto contFill = 0; // Fill style
+
+ // *********************************
+ // ** Hit occupancies
+
+ // ** Occupancy in XY plane
+ auto* pc_hit_ypos_vs_xpos = MakeCanvas<CbmQaCanvas>("hit_ypos_vs_xpos", "", 1600, 800);
+ pc_hit_ypos_vs_xpos->Divide2D(nSt);
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ pc_hit_ypos_vs_xpos->cd(iSt + 1);
+ fvph_hit_ypos_vs_xpos[iSt]->DrawCopy("colz", "");
+
+
+ // Reference size of the station saved to the detector interface
+ auto* pInnerCircle = new TEllipse(0., 0., fpDetInterface->GetRmin(iSt));
+ pInnerCircle->SetLineWidth(contWidth);
+ pInnerCircle->SetLineStyle(contStyle);
+ pInnerCircle->SetLineColor(contColor);
+ pInnerCircle->SetFillStyle(contFill);
+ pInnerCircle->Draw("SAME");
+
+ auto* pOuterCircle = new TEllipse(0., 0., fpDetInterface->GetRmax(iSt));
+ pOuterCircle->SetLineWidth(contWidth);
+ pOuterCircle->SetLineStyle(contStyle);
+ pOuterCircle->SetLineColor(contColor);
+ pOuterCircle->SetFillStyle(contFill);
+ pOuterCircle->Draw("SAME");
+
+ // Square boarders of station, which are used for the field approximation
+ double stXmin = +fpDetInterface->GetXmax(iSt);
+ double stXmax = -fpDetInterface->GetXmax(iSt);
+ double stYmin = -fpDetInterface->GetYmax(iSt);
+ double stYmax = +fpDetInterface->GetYmax(iSt);
+
+ auto* pBox = new TBox(stXmin, stYmin, stXmax, stYmax);
+ pBox->SetLineWidth(contWidth);
+ pBox->SetLineStyle(contStyle);
+ pBox->SetLineColor(contColor);
+ pBox->SetFillStyle(contFill);
+ pBox->Draw("SAME");
+ }
+
+ // ----- Occupancy projections
+ auto* pc_hit_xpos = MakeCanvas<CbmQaCanvas>("hit_xpos", "", 1400, 800);
+ pc_hit_xpos->Divide2D(nSt);
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ pc_hit_xpos->cd(iSt + 1);
+ auto* phProj = fvph_hit_ypos_vs_xpos[iSt]->ProjectionY();
+ phProj->DrawCopy();
+ }
+
+ auto* pc_hit_ypos = MakeCanvas<CbmQaCanvas>("hit_ypos", "", 1400, 800);
+ pc_hit_ypos->Divide2D(nSt);
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ pc_hit_ypos->cd(iSt + 1);
+ auto* phProj = fvph_hit_ypos_vs_xpos[iSt]->ProjectionX();
+ phProj->DrawCopy();
+ }
+
+
+ // ** Occupancy in XZ plane
+ MakeCanvas<CbmQaCanvas>("hit_xpos_vs_zpos", "", 600, 400);
+ fvph_hit_xpos_vs_zpos[nSt]->DrawCopy("colz", "");
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ // Station positions in detector IFS
+ double stZmin = fpDetInterface->GetZ(iSt) - 0.5 * fpDetInterface->GetThickness(iSt);
+ double stZmax = fpDetInterface->GetZ(iSt) + 0.5 * fpDetInterface->GetThickness(iSt);
+ double stHmin = -fpDetInterface->GetRmax(iSt);
+ double stHmax = +fpDetInterface->GetRmax(iSt);
+
+ auto* pBox = new TBox(stZmin, stHmin, stZmax, stHmax);
+ pBox->SetLineWidth(contWidth);
+ pBox->SetLineStyle(contStyle);
+ pBox->SetLineColor(contColor);
+ pBox->SetFillStyle(contFill);
+ pBox->Draw("SAME");
+ }
+
+ // ** Occupancy in YZ plane
+ MakeCanvas<CbmQaCanvas>("hit_ypos_vs_zpos", "", 600, 400);
+ fvph_hit_ypos_vs_zpos[nSt]->DrawCopy("colz", "");
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ // Station positions in detector IFS
+ double stZmin = fpDetInterface->GetZ(iSt) - 0.5 * fpDetInterface->GetThickness(iSt);
+ double stZmax = fpDetInterface->GetZ(iSt) + 0.5 * fpDetInterface->GetThickness(iSt);
+ double stHmin = -fpDetInterface->GetRmax(iSt);
+ double stHmax = +fpDetInterface->GetRmax(iSt);
+
+ auto* pBox = new TBox(stZmin, stHmin, stZmax, stHmax);
+ pBox->SetLineWidth(contWidth);
+ pBox->SetLineStyle(contStyle);
+ pBox->SetLineColor(contColor);
+ pBox->SetFillStyle(contFill);
+ pBox->Draw("SAME");
+ }
+
+ // ************
+ // ************
+
+ if (IsMCUsed()) {
+
+ // **********************
+ // ** Point occupancies
+
+ // ** Occupancy in XY plane
+ auto* pc_point_ypos_vs_xpos = MakeCanvas<CbmQaCanvas>("point_ypos_vs_xpos", "", 1600, 800);
+ pc_point_ypos_vs_xpos->Divide2D(nSt);
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ pc_point_ypos_vs_xpos->cd(iSt + 1);
+ fvph_point_ypos_vs_xpos[iSt]->DrawCopy("colz", "");
+
+ // Reference size of the station saved to the detector interface
+ auto* pInnerCircle = new TEllipse(0., 0., fpDetInterface->GetRmin(iSt));
+ pInnerCircle->SetLineWidth(contWidth);
+ pInnerCircle->SetLineStyle(contStyle);
+ pInnerCircle->SetLineColor(contColor);
+ pInnerCircle->SetFillStyle(contFill);
+ pInnerCircle->Draw("SAME");
+
+ auto* pOuterCircle = new TEllipse(0., 0., fpDetInterface->GetRmax(iSt));
+ pOuterCircle->SetLineWidth(contWidth);
+ pOuterCircle->SetLineStyle(contStyle);
+ pOuterCircle->SetLineColor(contColor);
+ pOuterCircle->SetFillStyle(contFill);
+ pOuterCircle->Draw("SAME");
+
+ // Square boarders of station, which are used for the field approximation
+ double stXmin = +fpDetInterface->GetXmax(iSt);
+ double stXmax = -fpDetInterface->GetXmax(iSt);
+ double stYmin = -fpDetInterface->GetYmax(iSt);
+ double stYmax = +fpDetInterface->GetYmax(iSt);
+
+ auto* pBox = new TBox(stXmin, stYmin, stXmax, stYmax);
+ pBox->SetLineWidth(contWidth);
+ pBox->SetLineStyle(contStyle);
+ pBox->SetLineColor(contColor);
+ pBox->SetFillStyle(contFill);
+ pBox->Draw("SAME");
+ }
+
+ // ** Occupancy in XZ plane
+ MakeCanvas<CbmQaCanvas>("point_xpos_vs_zpos", "", 600, 400);
+ fvph_point_xpos_vs_zpos[nSt]->SetStats(false);
+ fvph_point_xpos_vs_zpos[nSt]->DrawCopy("colz", "");
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ // Station positions in detector IFS
+ double stZmin = fpDetInterface->GetZ(iSt) - 0.5 * fpDetInterface->GetThickness(iSt);
+ double stZmax = fpDetInterface->GetZ(iSt) + 0.5 * fpDetInterface->GetThickness(iSt);
+ double stHmin = -fpDetInterface->GetRmax(iSt);
+ double stHmax = +fpDetInterface->GetRmax(iSt);
+
+ auto* pBox = new TBox(stZmin, stHmin, stZmax, stHmax);
+ pBox->SetLineWidth(contWidth);
+ pBox->SetLineStyle(contStyle);
+ pBox->SetLineColor(contColor);
+ pBox->SetFillStyle(contFill);
+ pBox->Draw("SAME");
+ }
+
+ // ** Occupancy in YZ plane
+ MakeCanvas<CbmQaCanvas>("point_ypos_vs_zpos", "", 600, 400);
+ fvph_point_ypos_vs_zpos[nSt]->SetStats(false);
+ fvph_point_ypos_vs_zpos[nSt]->DrawCopy("colz", "");
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ // Station positions in detector IFS
+ double stZmin = fpDetInterface->GetZ(iSt) - 0.5 * fpDetInterface->GetThickness(iSt);
+ double stZmax = fpDetInterface->GetZ(iSt) + 0.5 * fpDetInterface->GetThickness(iSt);
+ double stHmin = -fpDetInterface->GetRmax(iSt);
+ double stHmax = +fpDetInterface->GetRmax(iSt);
+
+ auto* pBox = new TBox(stZmin, stHmin, stZmax, stHmax);
+ pBox->SetLineWidth(contWidth);
+ pBox->SetLineStyle(contStyle);
+ pBox->SetLineColor(contColor);
+ pBox->SetFillStyle(contFill);
+ pBox->Draw("SAME");
+ }
+
+ // **************
+ // ** Residuals
+
+ // x-coordinate
+ auto* pc_res_x = MakeCanvas<CbmQaCanvas>("res_x", "Residuals for x coordinate", 1600, 800);
+ pc_res_x->Divide2D(nSt);
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ pc_res_x->cd(iSt + 1);
+ fvph_res_x[iSt]->DrawCopy("", "");
+ }
+
+ // y-coordinate
+ auto* pc_res_y = MakeCanvas<CbmQaCanvas>("res_y", "Residuals for y coordinate", 1600, 800);
+ pc_res_y->Divide2D(nSt);
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ pc_res_y->cd(iSt + 1);
+ fvph_res_y[iSt]->DrawCopy("", "");
+ }
+
+ // time
+ auto* pc_res_t = MakeCanvas<CbmQaCanvas>("res_t", "Residuals for time", 1600, 800);
+ pc_res_t->Divide2D(nSt);
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ pc_res_t->cd(iSt + 1);
+ fvph_res_t[iSt]->DrawCopy("", "");
+ }
+
+
+ // **********
+ // ** Pulls
+
+ // x-coordinate
+ auto* pc_pull_x = MakeCanvas<CbmQaCanvas>("pull_x", "Pulls for x coordinate", 1600, 800);
+ pc_pull_x->Divide2D(nSt);
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ pc_pull_x->cd(iSt + 1);
+ fvph_pull_x[iSt]->DrawCopy("", "");
+ }
+
+ // y-coordinate
+ auto* pc_pull_y = MakeCanvas<CbmQaCanvas>("pull_y", "Pulls for y coordinate", 1600, 800);
+ pc_pull_y->Divide2D(nSt);
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ pc_pull_y->cd(iSt + 1);
+ fvph_pull_y[iSt]->DrawCopy("", "");
+ }
+
+ // time
+ auto* pc_pull_t = MakeCanvas<CbmQaCanvas>("pull_t", "Pulls for time", 1600, 800);
+ pc_pull_t->Divide2D(nSt);
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ pc_pull_t->cd(iSt + 1);
+ fvph_pull_t[iSt]->DrawCopy("", "");
+ }
+
+
+ // ************************************
+ // ** Hit reconstruction efficiencies
+
+ auto* pc_reco_eff_vs_r =
+ MakeCanvas<CbmQaCanvas>("reco_eff_vs_r", "Hit efficiencies wrt distance from center", 1600, 800);
+ pc_reco_eff_vs_r->Divide2D(nSt);
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ pc_reco_eff_vs_r->cd(iSt + 1);
+ fvpe_reco_eff_vs_r[iSt]->Paint("AP");
+ auto* pGr = dynamic_cast<TGraphAsymmErrors*>(fvpe_reco_eff_vs_r[iSt]->GetPaintedGraph());
+ if (!pGr) {
+ LOG(error) << fName << ": unable to get painted graph from efficiency " << fvpe_reco_eff_vs_xy[iSt]->GetName();
+ continue;
+ }
+ pGr->DrawClone("AP");
+ auto* pFit = (TF1*) pGr->FindObject("pol0");
+ if (pFit) { pFit->Draw("SAME"); }
+ }
+
+ auto* pc_reco_eff_vs_xy = MakeCanvas<CbmQaCanvas>("reco_eff_vs_xy", "Hit efficiencies wrt x and y", 1600, 800);
+ pc_reco_eff_vs_xy->Divide2D(nSt);
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ pc_reco_eff_vs_xy->cd(iSt + 1);
+ fvpe_reco_eff_vs_xy[iSt]->Paint("colz");
+ auto* pH2 = dynamic_cast<TH2F*>(fvpe_reco_eff_vs_xy[iSt]->GetPaintedHistogram());
+ if (!pH2) {
+ LOG(error) << fName << ": unable to get painted histogram from efficiency "
+ << fvpe_reco_eff_vs_xy[iSt]->GetName();
+ continue;
+ }
+ pH2->DrawCopy("colz", "");
+ }
+ }
+
+ return kSUCCESS;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+InitStatus CbmCaInputQaTrd::InitHistograms()
+{
+ int nSt = fpDetInterface->GetNtrackingStations();
+
+ // ----- Histograms initialization
+ fvph_hit_ypos_vs_xpos.resize(nSt + 2, nullptr);
+ fvph_hit_ypos_vs_zpos.resize(nSt + 2, nullptr);
+ fvph_hit_xpos_vs_zpos.resize(nSt + 2, nullptr);
+
+ fvph_hit_station_delta_z.resize(nSt);
+
+ fvph_hit_dx.resize(nSt + 2, nullptr);
+ fvph_hit_dy.resize(nSt + 2, nullptr);
+ fvph_hit_dt.resize(nSt + 2, nullptr);
+
+ for (int iSt = 0; iSt < nSt + 2; ++iSt) {
+ TString nsuff = ""; // Histogram name suffix
+ TString tsuff = ""; // Histogram title suffix
+
+ // Select group: indexes 0 .. iSt - 1 stand for single TRD stations, index iSt == nSt is TRD-1D
+ // and iSt == nSt + 1 is TRD-2D
+ if (iSt < nSt) {
+ nsuff = Form("_st%d", iSt);
+ tsuff = Form(" in TRD station %d", iSt);
+ }
+ else if (iSt == nSt) {
+ nsuff = "_1D";
+ tsuff = " in TRD-1D";
+ }
+ else {
+ nsuff = "_2D";
+ tsuff = " in TRD-2D";
+ }
+
+ TString sN = "";
+ TString sT = "";
+
+ // Hit occupancy
+ sN = (TString) "hit_ypos_vs_xpos" + nsuff;
+ sT = (TString) "Hit occupancy in xy-Plane" + tsuff + ";x_{hit} [cm];y_{hit} [cm]";
+ fvph_hit_ypos_vs_xpos[iSt] = MakeHisto<TH2F>(sN, sT, kNbins, kRHitX[0], kRHitX[1], kNbins, kRHitY[0], kRHitY[1]);
+
+ sN = (TString) "hit_xpos_vs_zpos" + nsuff;
+ sT = (TString) "Hit occupancy in xz-Plane" + tsuff + ";z_{hit} [cm];x_{hit} [cm]";
+ fvph_hit_xpos_vs_zpos[iSt] = MakeHisto<TH2F>(sN, sT, kNbinsZ, kRHitZ[0], kRHitZ[1], kNbins, kRHitX[0], kRHitX[1]);
+
+ sN = (TString) "hit_ypos_vs_zpos" + nsuff;
+ sT = (TString) "Hit occupancy in yz-plane" + tsuff + ";z_{hit} [cm];y_{hit} [cm]";
+ fvph_hit_ypos_vs_zpos[iSt] = MakeHisto<TH2F>(sN, sT, kNbinsZ, kRHitZ[0], kRHitZ[1], kNbins, kRHitY[0], kRHitY[1]);
+
+ // Hit errors
+ sN = (TString) "hit_dx" + nsuff;
+ sT = (TString) "Hit position error along x-axis" + tsuff + ";dx_{hit} [cm]";
+ fvph_hit_dx[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, kRHitDx[0], kRHitDx[1]);
+
+ sN = (TString) "hit_dy" + nsuff;
+ sT = (TString) "Hit position error along y-axis" + tsuff + ";dy_{hit} [cm]";
+ fvph_hit_dy[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, kRHitDy[0], kRHitDy[1]);
+
+ sN = (TString) "hit_dt" + nsuff;
+ sT = (TString) "Hit time error" + tsuff + ";dt_{hit} [ns]";
+ fvph_hit_dt[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, kRHitDt[0], kRHitDt[1]);
+
+ if (iSt >= nSt) { continue; } // Following histograms are defined only for separate station
+
+ sN = (TString) "hit_station_delta_z" + nsuff;
+ sT = (TString) "Different between hit and station z-positions" + tsuff + ";z_{hit} - z_{st} [cm]";
+ fvph_hit_station_delta_z[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, -5., 5);
+
+ } // loop over station index: end
+
+ // ----- Initialize histograms, which are use MC-information
+ if (IsMCUsed()) {
+ // Resize histogram vectors
+ fvph_n_points_per_hit.resize(nSt + 2, nullptr);
+ fvph_point_ypos_vs_xpos.resize(nSt + 2, nullptr);
+ fvph_point_xpos_vs_zpos.resize(nSt + 2, nullptr);
+ fvph_point_ypos_vs_zpos.resize(nSt + 2, nullptr);
+ fvph_point_hit_delta_z.resize(nSt + 2, nullptr);
+ fvph_res_x.resize(nSt + 2, nullptr);
+ fvph_res_y.resize(nSt + 2, nullptr);
+ fvph_res_t.resize(nSt + 2, nullptr);
+ fvph_pull_x.resize(nSt + 2, nullptr);
+ fvph_pull_y.resize(nSt + 2, nullptr);
+ fvph_pull_t.resize(nSt + 2, nullptr);
+ fvph_res_x_vs_x.resize(nSt + 2, nullptr);
+ fvph_res_y_vs_y.resize(nSt + 2, nullptr);
+ fvph_res_t_vs_t.resize(nSt + 2, nullptr);
+ fvph_pull_x_vs_x.resize(nSt + 2, nullptr);
+ fvph_pull_y_vs_y.resize(nSt + 2, nullptr);
+ fvph_pull_t_vs_t.resize(nSt + 2, nullptr);
+ fvpe_reco_eff_vs_xy.resize(nSt + 2, nullptr);
+ fvpe_reco_eff_vs_r.resize(nSt + 2, nullptr);
+
+ for (int iSt = 0; iSt <= nSt; ++iSt) {
+ TString nsuff = ""; // Histogram name suffix
+ TString tsuff = ""; // Histogram title suffix
+
+ // Select group: indexes 0 .. iSt - 1 stand for single TRD stations, index iSt == nSt is TRD-1D
+ // and iSt == nSt + 1 is TRD-2D
+ if (iSt < nSt) {
+ nsuff = Form("_st%d", iSt);
+ tsuff = Form(" in TRD station %d", iSt);
+ }
+ else if (iSt == nSt) {
+ nsuff = "_1D";
+ tsuff = " in TRD-1D";
+ }
+ else {
+ nsuff = "_2D";
+ tsuff = " in TRD-2D";
+ }
+
+ TString sN = "";
+ TString sT = "";
+
+ sN = (TString) "n_points_per_hit" + nsuff;
+ sT = (TString) "Number of points per hit" + tsuff + ";N_{point}/hit";
+ fvph_n_points_per_hit[iSt] = MakeHisto<TH1F>(sN, sT, 10, -0.5, 9.5);
+
+ // Point occupancy
+ sN = (TString) "point_ypos_vs_xpos" + nsuff;
+ sT = (TString) "Point occupancy in XY plane" + tsuff + ";x_{MC} [cm];y_{MC} [cm]";
+ fvph_point_ypos_vs_xpos[iSt] =
+ MakeHisto<TH2F>(sN, sT, kNbins, kRHitX[0], kRHitX[1], kNbins, kRHitY[0], kRHitY[1]);
+
+ sN = (TString) "point_xpos_vs_zpos" + nsuff;
+ sT = (TString) "Point Occupancy in XZ plane" + tsuff + ";z_{MC} [cm];x_{MC} [cm]";
+ fvph_point_xpos_vs_zpos[iSt] =
+ MakeHisto<TH2F>(sN, sT, kNbinsZ, kRHitZ[0], kRHitZ[1], kNbins, kRHitX[0], kRHitX[1]);
+
+ sN = (TString) "point_ypos_vs_zpos" + nsuff;
+ sT = (TString) "Point Occupancy in YZ Plane" + tsuff + ";z_{MC} [cm];y_{MC} [cm]";
+ fvph_point_ypos_vs_zpos[iSt] =
+ MakeHisto<TH2F>(sN, sT, kNbinsZ, kRHitZ[0], kRHitZ[1], kNbins, kRHitY[0], kRHitY[1]);
+
+ // Difference between MC input z and hit z coordinates
+ sN = (TString) "point_hit_delta_z" + nsuff;
+ sT = (TString) "Distance between point and hit along z axis" + tsuff + ";z_{reco} - z_{MC} [cm]";
+ fvph_point_hit_delta_z[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, -0.04, 0.04);
+
+ sN = (TString) "res_x" + nsuff;
+ sT = (TString) "Residuals for x-coordinate" + tsuff + ";x_{reco} - x_{MC} [cm]";
+ fvph_res_x[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, kRResX[0], kRResX[1]);
+
+ sN = (TString) "res_y" + nsuff;
+ sT = (TString) "Residuals for y-coordinate" + tsuff + ";y_{reco} - y_{MC} [cm]";
+ fvph_res_y[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, kRResY[0], kRResY[1]);
+
+ sN = (TString) "res_t" + nsuff;
+ sT = (TString) "Residuals for time" + tsuff + ";t_{reco} - t_{MC} [ns]";
+ fvph_res_t[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, kRResT[0], kRResT[1]);
+
+ sN = (TString) "pull_x" + nsuff;
+ sT = (TString) "Pulls for x-coordinate" + tsuff + ";(x_{reco} - x_{MC}) / #sigma_{x}^{reco}";
+ fvph_pull_x[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, kRPullX[0], kRPullX[1]);
+
+ sN = (TString) "pull_y" + nsuff;
+ sT = (TString) "Pulls for y-coordinate" + tsuff + ";(y_{reco} - y_{MC}) / #sigma_{y}^{reco}";
+ fvph_pull_y[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, kRPullY[0], kRPullY[1]);
+
+ sN = (TString) "pull_t" + nsuff;
+ sT = (TString) "Pulls for time" + tsuff + ";(t_{reco} - t_{MC}) / #sigma_{t}^{reco}";
+ fvph_pull_t[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, kRPullT[0], kRPullT[1]);
+
+ sN = (TString) "res_x_vs_x" + nsuff;
+ sT = (TString) "Residuals for x-coordinate" + tsuff + ";x_{reco} [cm];x_{reco} - x_{MC} [cm]";
+ fvph_res_x_vs_x[iSt] = MakeHisto<TH2F>(sN, sT, kNbins, 0, 0, kNbins, kRResX[0], kRResX[1]);
+
+ sN = (TString) "res_y_vs_y" + nsuff;
+ sT = (TString) "Residuals for y-coordinate" + tsuff + ";y_{reco} [cm];y_{reco} - y_{MC} [cm]";
+ fvph_res_y_vs_y[iSt] = MakeHisto<TH2F>(sN, sT, kNbins, 0, 0, kNbins, kRResY[0], kRResY[1]);
+
+ sN = (TString) "res_t_vs_t" + nsuff;
+ sT = (TString) "Residuals for time" + tsuff + "t_{reco} [ns];t_{reco} - t_{MC} [ns]";
+ fvph_res_t_vs_t[iSt] = MakeHisto<TH2F>(sN, sT, kNbins, 0, 0, kNbins, kRResT[0], kRResT[1]);
+
+ sN = (TString) "pull_x_vs_x" + nsuff;
+ sT = (TString) "Pulls for x-coordinate" + tsuff + "x_{reco} [cm];(x_{reco} - x_{MC}) / #sigma_{x}^{reco}";
+ fvph_pull_x_vs_x[iSt] = MakeHisto<TH2F>(sN, sT, kNbins, 0, 0, kNbins, kRPullX[0], kRPullX[1]);
+
+ sN = (TString) "pull_y_vs_y" + nsuff;
+ sT = (TString) "Pulls for y-coordinate" + tsuff + ";y_{reco} [cm];(y_{reco} - y_{MC}) / #sigma_{y}^{reco}";
+ fvph_pull_y_vs_y[iSt] = MakeHisto<TH2F>(sN, sT, kNbins, 0, 0, kNbins, kRPullY[0], kRPullY[1]);
+
+ sN = (TString) "pull_t_vs_t" + nsuff;
+ sT = (TString) "Pulls for time" + tsuff + ";t_{reco} [ns];(t_{reco} - t_{MC}) / #sigma_{t}^{reco}";
+ fvph_pull_t_vs_t[iSt] = MakeHisto<TH2F>(sN, sT, kNbins, 0, 0, kNbins, kRPullT[0], kRPullT[1]);
+
+
+ sN = (TString) "reco_eff_vs_xy" + nsuff;
+ sT = (TString) "Hit rec. efficiency" + tsuff + ";x_{MC} [cm];y_{MC} [cm];#epsilon";
+ fvpe_reco_eff_vs_xy[iSt] = MakeEfficiency<CbmQaEff>(sN, sT, 100, -50, 50, 100, -50, 50);
+
+ sN = (TString) "reco_eff_vs_r" + nsuff;
+ sT = (TString) "Hit rec. efficiency" + tsuff + ";r_{MC} [cm];#epsilon";
+ fvpe_reco_eff_vs_r[iSt] = MakeEfficiency<CbmQaEff>(sN, sT, 100, 0, 100);
+ }
+ }
+ return kSUCCESS;
+}
diff --git a/reco/L1/qa/CbmCaInputQaTrd.h b/reco/L1/qa/CbmCaInputQaTrd.h
new file mode 100644
index 0000000000..9e9f403ea8
--- /dev/null
+++ b/reco/L1/qa/CbmCaInputQaTrd.h
@@ -0,0 +1,167 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// \file CbmCaInputQaTrd.h
+/// \date 13.01.2023
+/// \brief QA-task for CA tracking input from TRD detector (header)
+/// \author S.Zharko <s.zharko@gsi.de>
+
+
+#ifndef CbmCaInputQaTrd_h
+#define CbmCaInputQaTrd_h 1
+
+#include "CbmCaInputQaBase.h"
+#include "CbmMCDataManager.h"
+#include "CbmQaTask.h"
+
+#include "TMath.h"
+
+#include <set>
+#include <unordered_map>
+#include <vector>
+
+class CbmMCEventList;
+class CbmMCDataArray;
+class CbmMCDataManager;
+class CbmTimeSlice;
+class CbmMatch;
+class CbmTrdHit;
+class CbmTrdTrackingInterface;
+class TClonesArray;
+class TH1F;
+class TH2F;
+class CbmQaEff;
+
+/// A QA-task class, which provides assurance of TRD hits and geometry
+class CbmCaInputQaTrd : public CbmQaTask, public CbmCaInputQaBase {
+public:
+ /// Constructor from parameters
+ /// \param verbose Verbose level
+ /// \param isMCUsed Flag, whether MC information is available for this task
+ CbmCaInputQaTrd(int verbose, bool isMCUsed);
+
+ ClassDef(CbmCaInputQaTrd, 0);
+
+protected:
+ // ********************************************
+ // ** Virtual method override from CbmQaTask **
+ // ********************************************
+
+ /// Checks results of the QA
+ /// \return Success flag
+ bool Check();
+
+ /// Initializes data branches
+ InitStatus InitDataBranches();
+
+ /// Initializes canvases
+ InitStatus InitCanvases();
+
+ /// Initializes histograms
+ InitStatus InitHistograms();
+
+ /// Fills histograms per event or time-slice
+ void FillHistograms();
+
+ /// De-initializes histograms
+ void DeInit();
+
+private:
+ // *********************
+ // ** Private methods **
+ // *********************
+
+ // *********************
+ // ** Class variables **
+ // *********************
+
+
+ // ** Data branches **
+ CbmTrdTrackingInterface* fpDetInterface = nullptr; ///< Instance of detector interface
+
+ CbmTimeSlice* fpTimeSlice = nullptr; ///< Pointer to current time-slice
+
+ TClonesArray* fpHits = nullptr; ///< Array of hits
+
+ CbmMCDataManager* fpMCDataManager = nullptr; ///< MC data manager
+ CbmMCEventList* fpMCEventList = nullptr; ///< MC event list
+ CbmMCDataArray* fpMCTracks = nullptr; ///< Array of MC tracks
+ CbmMCDataArray* fpMCPoints = nullptr; ///< Array of MC points
+
+ TClonesArray* fpHitMatches = nullptr; ///< Array of hit matches
+
+ std::vector<char> fvTrdStationType; ///< Station type: TRD-1D and TRD-2D
+
+ // ** Histograms **
+ static constexpr double kEffRangeMin = 35.; ///< Lower limit of hit distance for efficiency integration [cm]
+ static constexpr double kEffRangeMax = 70.; ///< Upper limit of hit distance for efficiency integration [cm]
+
+ static constexpr double kRHitX[2] = {-100., 100}; ///< Range for hit x coordinate [cm]
+ static constexpr double kRHitY[2] = {-100., 100}; ///< Range for hit y coordinate [cm]
+ static constexpr double kRHitZ[2] = {30., 300.}; ///< Range for hit z coordinate [cm]
+
+ static constexpr int kNbins = 100; ///< General number of bins
+ static constexpr int kNbinsZ = 800; ///< Number of bins for z coordinate axis
+
+ static constexpr double kRHitDx[2] = {-.05, .005}; ///< Range for hit x coordinate [cm]
+ static constexpr double kRHitDy[2] = {-.05, .005}; ///< Range for hit y coordinate [cm]
+ static constexpr double kRHitDt[2] = {-10., 10.}; ///< Range for hit time [ns]
+
+ static constexpr double kRResX[2] = {-.4, .4}; ///< Range for residual of x coordinate [cm]
+ static constexpr double kRResY[2] = {-.4, .4}; ///< Range for residual of y coordinate [cm]
+ static constexpr double kRResT[2] = {-5., 5.}; ///< Range for residual of time [ns]
+
+ static constexpr double kRPullX[2] = {-10., 10.}; ///< Range for pull of x coordinate
+ static constexpr double kRPullY[2] = {-10., 10.}; ///< Range for pull of y coordinate
+ static constexpr double kRPullT[2] = {-5., 5.}; ///< Range for pull of time
+
+ // NOTE: the last three elements of each vector stands for integral distribution over TRD-1D and TRD-2D
+
+ // hit occupancy
+ std::vector<TH2F*> fvph_hit_ypos_vs_xpos; ///< hit ypos vs xpos in different stations
+ std::vector<TH2F*> fvph_hit_xpos_vs_zpos; ///< hit xpos vs zpos in different stations
+ std::vector<TH2F*> fvph_hit_ypos_vs_zpos; ///< hit ypos vs zpos in different stations
+
+ // difference between hit and station z positions
+ std::vector<TH1F*> fvph_hit_station_delta_z; ///< Difference between station and hit z positions [cm]
+
+ // hit errors
+ std::vector<TH1F*> fvph_hit_dx;
+ std::vector<TH1F*> fvph_hit_dy;
+ std::vector<TH1F*> fvph_hit_dt;
+
+ // MC points occupancy
+ std::vector<TH1F*> fvph_n_points_per_hit; ///< number of points per one hit in different stations
+
+
+ std::vector<TH2F*> fvph_point_ypos_vs_xpos; ///< point ypos vs xpos in different stations
+ std::vector<TH2F*> fvph_point_xpos_vs_zpos; ///< point xpos vs zpos in different stations
+ std::vector<TH2F*> fvph_point_ypos_vs_zpos; ///< point ypos vs zpos in different stations
+
+ std::vector<TH1F*> fvph_point_hit_delta_z; ///< difference between z of hit and MC point in different stations
+
+ // Residuals
+ std::vector<TH1F*> fvph_res_x; ///< residual for x coordinate in different stations
+ std::vector<TH1F*> fvph_res_y; ///< residual for y coordinate in different stations
+ std::vector<TH1F*> fvph_res_t; ///< residual for t coordinate in different stations
+
+ std::vector<TH2F*> fvph_res_x_vs_x; ///< residual for x coordinate in different stations
+ std::vector<TH2F*> fvph_res_y_vs_y; ///< residual for y coordinate in different stations
+ std::vector<TH2F*> fvph_res_t_vs_t; ///< residual for t coordinate in different stations
+
+ // Pulls
+ std::vector<TH1F*> fvph_pull_x; ///< pull for x coordinate in different stations
+ std::vector<TH1F*> fvph_pull_y; ///< pull for y coordinate in different stations
+ std::vector<TH1F*> fvph_pull_t; ///< pull for t coordinate in different stations
+
+ std::vector<TH2F*> fvph_pull_x_vs_x; ///< pull for x coordinate in different stations
+ std::vector<TH2F*> fvph_pull_y_vs_y; ///< pull for y coordinate in different stations
+ std::vector<TH2F*> fvph_pull_t_vs_t; ///< pull for t coordinate in different stations
+
+ // Hit efficiencies
+ std::vector<CbmQaEff*> fvpe_reco_eff_vs_xy; ///< Efficiency of hit reconstruction vs. x and y coordinates of MC
+ std::vector<CbmQaEff*> fvpe_reco_eff_vs_r; ///< Efficiency of hit reconstruction vs. distance from center
+};
+
+#endif // CbmCaInputQaTrd_h
diff --git a/reco/L1/qa/CbmCaQaCuts.h b/reco/L1/qa/CbmCaQaCuts.h
index 08810ac00a..9e9d643ec7 100644
--- a/reco/L1/qa/CbmCaQaCuts.h
+++ b/reco/L1/qa/CbmCaQaCuts.h
@@ -7,8 +7,7 @@
namespace cbm::ca::qa::cuts
{
- constexpr double kMinP = 0.05; ///< minimal momentum [Gev/c]
- constexpr double kMinEff = 0.50; ///< minimal efficiency of hit finders (NOTE: value is tmp)
+ constexpr double kMinP = 0.05; ///< minimal momentum [Gev/c]
// Max difference between hit and station z
constexpr double kMaxDzStHitSts = 1.; ///< Max distance for STS
diff --git a/reco/L1/qa/CbmTofInteraction.cxx b/reco/L1/qa/CbmTofInteraction.cxx
new file mode 100644
index 0000000000..f6ecac3db7
--- /dev/null
+++ b/reco/L1/qa/CbmTofInteraction.cxx
@@ -0,0 +1,99 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// \file CbmTofInteraction.h
+/// \date 02.02.2023
+/// \brief Representation of MC track interaction with a TOF module
+/// \author P.-A. Loizeau
+/// \author S. Zharko
+
+
+#include "CbmTofInteraction.h"
+
+#include "Logger.h"
+
+#include <sstream>
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+CbmTofInteraction::CbmTofInteraction() : CbmTofPoint(), fNofPoints(0) {};
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void CbmTofInteraction::AddPoint(const CbmTofPoint* pPoint)
+{
+ // ----- Update track and detector ID
+ if (fNofPoints == 0) {
+ fTrackID = pPoint->GetTrackID();
+ fDetectorID = pPoint->GetDetectorID();
+ }
+ else {
+ // Track and detector IDs schold be equal for every point
+ LOG_IF(fatal, pPoint->GetTrackID() != fTrackID) << "Attempt to add point with inconsistent track ID";
+ LOG_IF(fatal, pPoint->GetDetectorID() != fDetectorID) << "Attempt to add point with inconsistent detector ID";
+ }
+
+ // ----- Update position and momenta
+ UpdateAverage(pPoint->GetX(), fX);
+ UpdateAverage(pPoint->GetY(), fY);
+ UpdateAverage(pPoint->GetZ(), fZ);
+ UpdateAverage(pPoint->GetTime(), fTime);
+
+ UpdateAverage(pPoint->GetPx(), fPx);
+ UpdateAverage(pPoint->GetPy(), fPy);
+ UpdateAverage(pPoint->GetPz(), fPz);
+ UpdateAverage(pPoint->GetLength(), fLength); // Is it correct?
+ fELoss += pPoint->GetEnergyLoss(); // Is it correct?
+ fvpPoints.push_back(pPoint);
+ fNofPoints++;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void CbmTofInteraction::Clear()
+{
+ fTrackID = -1;
+ fDetectorID = -1;
+ fEventId = 0;
+ fPx = 0.;
+ fPy = 0.;
+ fPz = 0.;
+ fTime = 0.;
+ fLength = 0.;
+ fELoss = 0.;
+ fX = 0.;
+ fY = 0.;
+ fZ = 0.;
+ fNofPoints = 0;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void CbmTofInteraction::SetFromPoint(const CbmTofPoint* pPoint)
+{
+ // Check consistency of the point and the interaction
+ LOG_IF(fatal, fTrackID != pPoint->GetTrackID() || fDetectorID != pPoint->GetDetectorID())
+ << "CbmTofInteraction: attempt to add point with inconsistent track or detector IDs: track " << fTrackID << " vs. "
+ << pPoint->GetTrackID() << ", sensor " << fDetectorID << " vs. " << pPoint->GetDetectorID();
+
+ fX = pPoint->GetX();
+ fY = pPoint->GetY();
+ fZ = pPoint->GetZ();
+ fTime = pPoint->GetTime();
+ fPx = pPoint->GetPx();
+ fPy = pPoint->GetPy();
+ fPz = pPoint->GetPz();
+ fLength = pPoint->GetLength();
+ fELoss = pPoint->GetEnergyLoss();
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+std::string CbmTofInteraction::ToString() const
+{
+ std::stringstream msg;
+ msg << CbmTofPoint::ToString();
+ msg << " Number of points: " << fNofPoints;
+ return msg.str();
+}
diff --git a/reco/L1/qa/CbmTofInteraction.h b/reco/L1/qa/CbmTofInteraction.h
new file mode 100644
index 0000000000..33a65b8c6e
--- /dev/null
+++ b/reco/L1/qa/CbmTofInteraction.h
@@ -0,0 +1,77 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// \file CbmTofInteraction.h
+/// \date 02.02.2023
+/// \brief Representation of MC track interaction with a TOF module
+/// \author P.-A. Loizeau
+/// \author S. Zharko
+
+#ifndef CbmTofInteraction_h
+#define CbmTofInteraction_h 1
+
+#include "CbmTofPoint.h"
+
+#include <string>
+#include <vector>
+/// Class describes an interaction of a MC track with a TOF module as a solid structure
+///
+class CbmTofInteraction : public CbmTofPoint {
+public:
+ // ----- Constructors and destructor
+ /// Default constructor
+ CbmTofInteraction();
+
+ /// Constructor with signature from CbmTofPoint
+ /// \param args List of parameters for any CbmTofPoint constructor
+ template<typename... Args>
+ CbmTofInteraction(Args... args) : CbmTofPoint(args...)
+ , fNofPoints(0)
+ {
+ }
+
+ /// Destructor
+ ~CbmTofInteraction() = default;
+
+ // ----- Modifiers
+ /// \brief Adds a point to the interaction
+ /// New point updates the following properties of the interaction: position, entrance momentum
+ /// \param pPoint Pointer to TOF MC-point
+ void AddPoint(const CbmTofPoint* pPoint);
+
+ /// \brief Clears the instance
+ void Clear();
+
+ /// \brief Sets parameters from a particular TOF MC-point
+ /// \param pPoint Pointer to TOF point
+ /// The purpose of the
+ void SetFromPoint(const CbmTofPoint* pPoint);
+
+ // ----- Getters
+ /// \brief Gets number of stored points
+ int GetNofPoints() const { return fNofPoints; }
+
+ /// \brief Saves content of the class to string
+ std::string ToString() const;
+
+ /// \brief Gets pointers to points (TMP!!!!)
+ const std::vector<const CbmTofPoint*>& GetPoints() const { return fvpPoints; }
+
+private:
+ // ----- Utility functions
+ /// \brief Updates average of the property from original TOF MC-point
+ /// \param update Property of the added MC point
+ /// \param property Updated property of this interaction
+ template<typename T>
+ void UpdateAverage(const T& update, T& property)
+ {
+ property = (fNofPoints * property + update) / (fNofPoints + 1);
+ }
+
+ // ----- Properties
+ int fNofPoints = 0; ///< Number of CbmTofPoint objects, from which the interaction is constructed
+ std::vector<const CbmTofPoint*> fvpPoints; ///< Vector of point pointer (TMP!!!!)
+};
+
+#endif // CbmTofInteraction_h
--
GitLab