From 302977ed9bc2ca73c4642c9b7272fd7766977dd4 Mon Sep 17 00:00:00 2001
From: "s.zharko@gsi.de" <s.zharko@gsi.de>
Date: Mon, 23 Jan 2023 17:37:19 +0100
Subject: [PATCH] CA: QA task for tracking input from STS (v.2)
---
core/qa/CMakeLists.txt | 6 +
core/qa/CbmQaBaseLinkDef.h | 2 +
core/qa/CbmQaCanvas.h | 11 +
core/qa/CbmQaHist.h | 19 +-
core/qa/CbmQaTask.cxx | 170 ++++
core/qa/CbmQaTask.h | 229 +++++
macro/run/run_reco_qa.C | 292 +++++++
reco/L1/CMakeLists.txt | 3 +
reco/L1/CbmCaMCModule.h | 507 +++++++++++
reco/L1/L1LinkDef.h | 2 +
reco/L1/qa/CbmCaInputQaSts.cxx | 1193 ++++++++++++++++++++++++++
reco/L1/qa/CbmCaInputQaSts.h | 201 +++++
reco/L1/qa/CbmCaQaCuts.h | 17 +
reco/L1/qa/CbmTrackerInputQaTrd.cxx | 2 +
reco/L1/qa/CbmTrackingInputQaSts.cxx | 20 +-
15 files changed, 2643 insertions(+), 31 deletions(-)
create mode 100644 core/qa/CbmQaTask.cxx
create mode 100644 core/qa/CbmQaTask.h
create mode 100644 macro/run/run_reco_qa.C
create mode 100644 reco/L1/CbmCaMCModule.h
create mode 100644 reco/L1/qa/CbmCaInputQaSts.cxx
create mode 100644 reco/L1/qa/CbmCaInputQaSts.h
create mode 100644 reco/L1/qa/CbmCaQaCuts.h
diff --git a/core/qa/CMakeLists.txt b/core/qa/CMakeLists.txt
index 5a87997b4e..3589694a85 100644
--- a/core/qa/CMakeLists.txt
+++ b/core/qa/CMakeLists.txt
@@ -7,12 +7,14 @@ set(SRCS
CbmQaPie.cxx
CbmQaHist.cxx
CbmQaTable.cxx
+ CbmQaTask.cxx
)
set(LIBRARY_NAME CbmQaBase)
set(LINKDEF ${LIBRARY_NAME}LinkDef.h)
set(PUBLIC_DEPENDENCIES
+ FairRoot::Base
FairLogger::FairLogger
ROOT::Core
ROOT::Gpad
@@ -24,3 +26,7 @@ set(INTERFACE_DEPENDENCIES
)
generate_cbm_library()
+
+Install(FILES CbmQaTask.h CbmQaCanvas.h CbmQaTable.h CbmQaHist.h
+ DESTINATION include
+ )
diff --git a/core/qa/CbmQaBaseLinkDef.h b/core/qa/CbmQaBaseLinkDef.h
index 9da99d54d5..5d2349aa39 100644
--- a/core/qa/CbmQaBaseLinkDef.h
+++ b/core/qa/CbmQaBaseLinkDef.h
@@ -16,6 +16,7 @@
#pragma link C++ class CbmQaPie - ;
// create streamers automatically
+#pragma link C++ class CbmQaEff + ;
#pragma link C++ class CbmQaPieSlice + ;
#pragma link C++ class CbmQaHist < TH1F> + ;
#pragma link C++ class CbmQaHist < TH1D> + ;
@@ -23,5 +24,6 @@
#pragma link C++ class CbmQaHist < TProfile> + ;
#pragma link C++ class CbmQaHist < TProfile2D> + ;
#pragma link C++ class CbmQaTable + ;
+#pragma link C++ class CbmQaTask + ;
#endif
diff --git a/core/qa/CbmQaCanvas.h b/core/qa/CbmQaCanvas.h
index ffcf663812..bce7112517 100644
--- a/core/qa/CbmQaCanvas.h
+++ b/core/qa/CbmQaCanvas.h
@@ -46,6 +46,17 @@ public:
/// Divide canvas into nPads in 2D in a nice way
void Divide2D(int nPads);
+
+ /// a static canvas for temporary drawing
+ static CbmQaCanvas& GetDummyCanvas()
+ {
+ /// the static variable will be initialised at the first call;
+ /// deleted at the application end (c++11)
+ static CbmQaCanvas dummy("CbmQaTempCanvas", "CbmQaTempCanvas", 1, 1);
+ return dummy;
+ }
+
+
private:
/// Use a specific type name in order to avoid ambiguities when unrolling
/// templates
diff --git a/core/qa/CbmQaHist.h b/core/qa/CbmQaHist.h
index 546e573b25..194710611c 100644
--- a/core/qa/CbmQaHist.h
+++ b/core/qa/CbmQaHist.h
@@ -70,10 +70,10 @@ public:
TFitResultPtr Fit(Types... args)
{
TVirtualPad* padsav = gPad;
- GetDummyCanvas().cd();
+ CbmQaCanvas::GetDummyCanvas().cd();
this->Sumw2();
auto ret = HistTypeT::Fit(args...);
- GetDummyCanvas().Clear();
+ CbmQaCanvas::GetDummyCanvas().Clear();
// make the output look nice
@@ -118,15 +118,15 @@ public:
int saveStat = gStyle->GetOptStat();
int saveFit = gStyle->GetOptFit();
- GetDummyCanvas().cd();
+ CbmQaCanvas::GetDummyCanvas().cd();
gStyle->SetOptStat(fOptStat);
gStyle->SetOptFit(fOptFit);
this->SetStats(0); // remove the old stat window
this->SetStats(1); // set the flag to create the stat window during Draw()
this->Draw();
- GetDummyCanvas().Update();
- GetDummyCanvas().Clear();
+ CbmQaCanvas::GetDummyCanvas().Update();
+ CbmQaCanvas::GetDummyCanvas().Clear();
// restore the environment
gStyle->SetOptStat(saveStat);
@@ -135,15 +135,6 @@ public:
}
private:
- /// a static canvas for temporary drawing
- static CbmQaCanvas& GetDummyCanvas()
- {
- /// the static variable will be initialised at the first call;
- /// deleted at the application end (c++11)
- static CbmQaCanvas dummy("CbmQaTempCanvas", "CbmQaTempCanvas", 1, 1);
- return dummy;
- }
-
int fOptStat = 1;
int fOptFit = 0;
diff --git a/core/qa/CbmQaTask.cxx b/core/qa/CbmQaTask.cxx
new file mode 100644
index 0000000000..d74cd4582c
--- /dev/null
+++ b/core/qa/CbmQaTask.cxx
@@ -0,0 +1,170 @@
+/* opyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// \file CbmQaTask.cxx
+/// \brief A base class for CBM QA task logic (implementation)
+/// \author S.Zharko <s.zharko@lx-pool.gsi.de>
+/// \data 12.01.2023
+
+
+#include "CbmQaTask.h"
+
+#include "CbmQaCanvas.h"
+
+#include "FairRootManager.h"
+#include "FairRunAna.h"
+#include "FairSink.h"
+
+#include "TAxis.h"
+
+#include <array>
+
+ClassImp(CbmQaTask);
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+CbmQaTask::CbmQaTask(const char* name, int verbose, bool isMCUsed) : FairTask(name, verbose), fbUseMC(isMCUsed) {}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void CbmQaTask::Exec(Option_t* /*option*/)
+{
+ fNofEvents.SetVal(fNofEvents.GetVal() + 1);
+ LOG_IF(info, fVerbose > 1) << fName << ": event " << fNofEvents.GetVal();
+ FillHistograms();
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void CbmQaTask::Finish()
+{
+ // Processes histograms in the end of the run
+ AnalyzeHistograms();
+
+ // Check QA
+ Check();
+
+ // Processes canvases in the end of the run
+ LOG_IF(info, fVerbose > 1) << fName << ": initializing canvases";
+ InitCanvases();
+
+ // 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);
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+InitStatus CbmQaTask::Init()
+{
+ LOG_IF(info, fVerbose > 0) << fName << ": initializing task ...";
+
+ InitStatus res = kSUCCESS;
+
+ // ----- Clear map of the histograms (note)
+ DeInitBase();
+
+ // ----- Initialize data branches
+ LOG_IF(info, fVerbose > 1) << fName << ": initializing input data branches";
+ res = std::max(res, InitDataBranches());
+
+ // ----- Initialize I/O
+ fpFolderRoot = std::make_unique<TFolder>(fName, fName); // The name of the folder follows the name of the task
+ fpFolderRoot->SetOwner(true); // When true, TFolder owns all added objects
+ fpFolderRoot->Add(&fNofEvents);
+
+ // ----- Initialize histograms and canvases
+ LOG_IF(info, fVerbose > 1) << fName << ": initializing histograms";
+ res = std::max(res, InitHistograms());
+
+ fNofEvents.SetVal(0);
+
+ return res;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+InitStatus CbmQaTask::ReInit()
+{
+ LOG_IF(info, fVerbose > 2) << fName << "::DeInitBase() is called";
+ return Init();
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void CbmQaTask::DeInitBase()
+{
+ LOG_IF(info, fVerbose > 2) << fName << "::DeInitBase() is called";
+ // De-initialize basic data members
+ // ...
+
+ delete fpFolderRoot.get();
+ fpFolderRoot = nullptr;
+
+ // De-initialize particular QA-task implementation
+ DeInit();
+}
+
+
+// ---------------------------------------------------------------------------------------------------------------------
+// NOTE: Example
+void CbmQaTask::SetHistoProperties(TH1* pHisto)
+{
+ constexpr double kHdefTextSize = 0.04;
+
+ pHisto->SetStats(kTRUE);
+ pHisto->Sumw2();
+
+ // Axis property settings
+ std::array<TAxis*, 3> vpAxes = {pHisto->GetXaxis(), pHisto->GetYaxis(), pHisto->GetZaxis()};
+ for (auto* pAxis : vpAxes) {
+ pAxis->SetTitleSize(kHdefTextSize);
+ pAxis->SetLabelSize(kHdefTextSize);
+ }
+ vpAxes[0]->SetNdivisions(504);
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+InitStatus CbmQaTask::InitDataBranches()
+{
+ LOG_IF(info, fVerbose > 1) << fName << ": data branches initialization function is not defined";
+ return kSUCCESS;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+InitStatus CbmQaTask::InitHistograms()
+{
+ LOG_IF(info, fVerbose > 1) << fName << ": histogram initialization function is not defined";
+ return kSUCCESS;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void CbmQaTask::FillHistograms()
+{
+ LOG_IF(info, fVerbose > 1) << fName << ": histogram filling function is not defined";
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void CbmQaTask::AnalyzeHistograms()
+{
+ LOG_IF(info, fVerbose > 1) << fName << ": no action on histograms in the end of the run is provided";
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+InitStatus CbmQaTask::InitCanvases()
+{
+ LOG_IF(info, fVerbose > 1) << fName << ": no initialization of canvases is provided";
+ return kSUCCESS;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void CbmQaTask::DeInit() { LOG_IF(info, fVerbose > 1) << fName << ": no extra de-initialization is provided"; }
diff --git a/core/qa/CbmQaTask.h b/core/qa/CbmQaTask.h
new file mode 100644
index 0000000000..c18e596485
--- /dev/null
+++ b/core/qa/CbmQaTask.h
@@ -0,0 +1,229 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// \file CbmQaTask.h
+/// \brief A base class for CBM QA task logic
+/// \author S.Zharko <s.zharko@gsi.de>
+/// \data 12.01.2022
+
+#ifndef CbmQaTask_h
+#define CbmQaTask_h 1
+
+#include "FairTask.h"
+#include "Logger.h"
+
+#include "TCanvas.h"
+#include "TEfficiency.h"
+#include "TFolder.h"
+#include "TH1.h"
+#include "TParameter.h"
+#include "TROOT.h"
+
+#include <string_view>
+
+/// Class CbmQaTask is to be inherited with a particular QA-task. It provides mechanisms for storage and management
+/// of QA canvases and histograms management
+class CbmQaTask : public FairTask {
+public:
+ /// Constructor from parameters
+ /// \param name Name of the task
+ /// \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);
+
+ /// Default constructor
+ CbmQaTask() = delete; // TODO: Let's see, what can happen, if one deletes default constructor
+
+ /// Destructor
+ virtual ~CbmQaTask() = default;
+
+ // Copy and move semantics
+ CbmQaTask(const CbmQaTask&) = delete;
+ CbmQaTask(CbmQaTask&&) = delete;
+ CbmQaTask& operator=(const CbmQaTask&) = delete;
+ CbmQaTask& operator=(CbmQaTask&&) = delete;
+
+ /// Returns flag, whether MC information is used or not in the task
+ bool IsMCUsed() const { return fbUseMC; }
+
+ /// FairTask: Task initialization in the beginning of the run
+ InitStatus Init();
+
+ /// FairTask: Task reinitialization
+ InitStatus ReInit();
+
+ /// FairTask: Defines action of the task in the event/TS
+ void Exec(Option_t* /*option*/);
+
+ /// FairTask: Defines action of the task in the end of run
+ void Finish();
+
+ ClassDef(CbmQaTask, 0);
+
+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; }
+
+ /// De-initialize the task
+ virtual void DeInit();
+
+ /// Initializes data branches
+ virtual InitStatus InitDataBranches();
+
+ /// Initializes histograms
+ virtual InitStatus InitHistograms();
+
+ /// Initializes canvases
+ virtual InitStatus InitCanvases();
+
+ /// 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
+ /// \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
+ /// \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
+ /// \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);
+
+ /// Method to setup histogram properties (text sizes etc.)
+ /// \param pHist Pointer to a histogram to tune
+ virtual void SetHistoProperties(TH1* pHist);
+
+ /// Gets a reference to histograms map
+ ///const auto& GetHistoMap() const { return fmpHistList; }
+
+ /// Gets a reference to canvases map
+ ///const auto& GetCanvasMap() const { return fmpCanvList; }
+
+ /// Get current event number
+ int GetEventNumber() const { return fNofEvents.GetVal(); }
+
+
+private:
+ // **********************************************************
+ // ** Functions accessible inside the CbmQaTask class only **
+ // **********************************************************
+
+ /// De-initialize this task
+ void DeInitBase();
+
+ // ***********************
+ // ** 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
+
+ TParameter<int> fNofEvents {"nEvents", 0}; ///< Number of processed events
+};
+
+
+// *************************************************
+// ** Inline and template function implementation **
+// *************************************************
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+template<typename T, typename... Args>
+T* CbmQaTask::MakeCanvas(const char* name, Args... args)
+{
+ if (gROOT->FindObject(name)) {
+ LOG(warn) << fName << ": A canvas with name \"" << name << "\" was previously created and will be deleted now "
+ << "to avoid memory leaks";
+ T* pCanv = (T*) gROOT->FindObject(name);
+ delete pCanv;
+ }
+
+ // Create a new canvas
+ T* pCanv = new T(name, args...);
+ pCanv->SetLeftMargin(0.2);
+ pCanv->SetBottomMargin(0.2);
+
+
+ // Register canvas in the folder
+ if (!fpFolderCanv) { fpFolderCanv = fpFolderRoot->AddFolder("canvases", "Canvases"); }
+ fpFolderCanv->Add(pCanv);
+
+ return pCanv;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+template<typename T, typename... Args>
+T* CbmQaTask::MakeEfficiency(const char* name, Args... args)
+{
+ if (gROOT->FindObject(name)) {
+ LOG(warn) << fName << ": An efficiency with name \"" << name << "\" was previously created and will be deleted now "
+ << "to avoid memory leaks";
+ auto* pEff = (T*) gROOT->FindObject(name);
+ delete pEff;
+ }
+
+ // Create a new efficiency
+ auto* pEff = new T(name, args...);
+
+ // Register efficiency in the folder
+ if (!fpFolderEff) { fpFolderEff = fpFolderRoot->AddFolder("efficiencies", "Efficiencies"); }
+ fpFolderEff->Add(pEff);
+
+ return pEff;
+}
+
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+template<typename T, typename... Args>
+T* CbmQaTask::MakeHisto(const char* name, Args... args)
+{
+ // Check, if the histogram with a given name was already created. If so, delete it
+ if (gROOT->FindObject(name)) {
+ LOG(warn) << fName << ": A histogram with name \"" << name << "\" was previously created and will be deleted now "
+ << "to avoid memory leaks";
+ T* pHist = (T*) gROOT->FindObject(name);
+ delete pHist;
+ }
+
+ // 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"); }
+ fpFolderHist->Add(pHist);
+
+ SetHistoProperties(pHist);
+
+ return pHist;
+}
+
+#endif // CbmQaTask_h
diff --git a/macro/run/run_reco_qa.C b/macro/run/run_reco_qa.C
new file mode 100644
index 0000000000..f286a031ea
--- /dev/null
+++ b/macro/run/run_reco_qa.C
@@ -0,0 +1,292 @@
+/* Copyright (C) 2006-2020 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergey Gorbunov, Denis Bertini [committer], Florian Uhlig */
+
+// --------------------------------------------------------------------------
+//
+// Macro for simulation & reconstruction QA
+//
+// The following naming conventions are assumed:
+// Raw data file: [dataset].event.raw.root
+// Transport file: [dataset].tra.root
+// Parameter file: [dataset].par.root
+// Reconstruction file: [dataset].rec.root
+//
+// S. Gorbunov 28/09/2020
+//
+// --------------------------------------------------------------------------
+
+// Includes needed for IDE
+#if !defined(__CLING__)
+#include "CbmDefs.h"
+#include "CbmMCDataManager.h"
+#include "CbmMuchDigitizerQa.h"
+#include "CbmMuchHitFinderQa.h"
+#include "CbmMuchTransportQa.h"
+#include "CbmSetup.h"
+
+#include <FairFileSource.h>
+#include <FairMonitor.h>
+#include <FairParAsciiFileIo.h>
+#include <FairParRootFileIo.h>
+#include <FairRootFileSink.h>
+#include <FairRunAna.h>
+#include <FairRuntimeDb.h>
+#include <FairSystemInfo.h>
+
+#include <TStopwatch.h>
+#endif
+
+void run_reco_qa(TString dataTra = "data/sis100_muon_jpsi_test", TString dataRaw = "data/sis100_muon_jpsi_test",
+ TString dataReco = "data/sis100_muon_jpsi_test", TString dataPar = "data/sis100_muon_jpsi_test",
+ TString dataSink = "data/sis100_muon_jpsi_test", TString setup = "sis100_muon_jpsi",
+ Int_t nEvents = -1, TString dataTra2 = "", TString dataTra3 = "")
+{
+
+ gROOT->SetBatch(kTRUE);
+
+ // ========================================================================
+ // Adjust this part according to your requirements
+
+ // ----- Logger settings ----------------------------------------------
+ FairLogger::GetLogger()->SetLogScreenLevel("INFO");
+ fair::Logger::DefineVerbosity(
+ "user1", fair::VerbositySpec::Make(fair::VerbositySpec::Info::severity, fair::VerbositySpec::Info::file_line));
+ FairLogger::GetLogger()->SetLogVerbosityLevel("user1");
+ FairLogger::GetLogger()->SetColoredLog(true);
+ // ------------------------------------------------------------------------
+
+ // ----- Environment --------------------------------------------------
+ int verbose = 100;
+ bool bUseMC = true;
+ TString myName = "run_qa"; // this macro's name for screen output
+ TString srcDir = gSystem->Getenv("VMCWORKDIR"); // top source directory
+ // ------------------------------------------------------------------------
+
+ // ----- In- and output file names ------------------------------------
+ TString traFile = dataTra + ".tra.root";
+ TString tra2File = dataTra2 + ".tra.root";
+ TString tra3File = dataTra3 + ".tra.root";
+ TString rawFile = dataRaw + ".raw.root";
+ TString parFile = dataPar + ".par.root";
+ TString recFile = dataReco + ".reco.root";
+ TString sinkFile = dataSink + ".qa.root";
+ // ------------------------------------------------------------------------
+
+ // ----- Load the geometry setup -------------------------------------
+ std::cout << std::endl;
+ std::cout << "-I- " << myName << ": Loading setup " << setup << std::endl;
+ CbmSetup::Instance()->LoadSetup(setup);
+ // You can modify the pre-defined setup by using
+ // CbmSetup::Instance()->RemoveModule(ESystemId) or
+ // CbmSetup::Instance()->SetModule(ESystemId, const char*, Bool_t) or
+ // CbmSetup::Instance()->SetActive(ESystemId, Bool_t)
+ // See the class documentation of CbmSetup.
+ // ------------------------------------------------------------------------
+
+ // ----- Parameter files as input to the runtime database -------------
+ std::cout << std::endl;
+ std::cout << "-I- " << myName << ": Defining parameter files " << std::endl;
+ TList* parFileList = new TList();
+ TString geoTag;
+
+ // - MUCH digitisation parameters
+ TString muchParFile {};
+ if (CbmSetup::Instance()->GetGeoTag(ECbmModuleId::kMuch, geoTag)) {
+ bool mcbmFlag = geoTag.Contains("mcbm", TString::kIgnoreCase);
+ muchParFile = srcDir + "/parameters/much/much_";
+ muchParFile += (mcbmFlag) ? geoTag : geoTag(0, 4);
+ muchParFile += "_digi_sector.root";
+ { // init geometry from the file
+ TFile* f = new TFile(muchParFile, "R");
+ TObjArray* stations = (TObjArray*) f->Get("stations");
+ CbmMuchGeoScheme::Instance()->Init(stations, mcbmFlag);
+ }
+ }
+
+ // - TRD digitisation parameters
+ if (CbmSetup::Instance()->GetGeoTag(ECbmModuleId::kTrd, geoTag)) {
+ const Char_t* npar[4] = {"asic", "digi", "gas", "gain"};
+ TObjString* trdParFile(NULL);
+ for (Int_t i(0); i < 4; i++) {
+ trdParFile = new TObjString(srcDir + "/parameters/trd/trd_" + geoTag + "." + npar[i] + ".par");
+ parFileList->Add(trdParFile);
+ std::cout << "-I- " << myName << ": Using parameter file " << trdParFile->GetString() << std::endl;
+ }
+ }
+
+ // - TOF digitisation parameters
+ if (CbmSetup::Instance()->GetGeoTag(ECbmModuleId::kTof, geoTag)) {
+ TObjString* tofBdfFile = new TObjString(srcDir + "/parameters/tof/tof_" + geoTag + ".digibdf.par");
+ parFileList->Add(tofBdfFile);
+ std::cout << "-I- " << myName << ": Using parameter file " << tofBdfFile->GetString() << std::endl;
+ }
+ // ------------------------------------------------------------------------
+
+ // In general, the following parts need not be touched
+ // ========================================================================
+
+ // ----- Timer --------------------------------------------------------
+ TStopwatch timer;
+ timer.Start();
+ // ------------------------------------------------------------------------
+
+ // ---- Debug option -------------------------------------------------
+ gDebug = 0;
+ // ------------------------------------------------------------------------
+
+ // ----- FairRunAna ---------------------------------------------------
+ FairFileSource* inputSource = new FairFileSource(rawFile);
+ inputSource->AddFriend(traFile);
+ inputSource->AddFriend(recFile);
+ if (!dataTra2.IsNull()) inputSource->AddFriend(tra2File);
+ if (!dataTra3.IsNull()) inputSource->AddFriend(tra3File);
+
+ FairRunAna* run = new FairRunAna();
+ run->SetSource(inputSource);
+ run->SetGenerateRunInfo(kFALSE);
+
+ FairRootFileSink* sink = new FairRootFileSink(sinkFile);
+ run->SetSink(sink);
+
+ TString monitorFile {sinkFile};
+ monitorFile.ReplaceAll("qa", "qa.monitor");
+ FairMonitor::GetMonitor()->EnableMonitor(kTRUE, monitorFile);
+ // ------------------------------------------------------------------------
+
+ // ----- MCDataManager -----------------------------------
+ if (bUseMC) {
+ CbmMCDataManager* mcManager = new CbmMCDataManager("MCDataManager", 0);
+ mcManager->AddFile(traFile);
+ if (!dataTra2.IsNull()) mcManager->AddFile(tra2File);
+ if (!dataTra3.IsNull()) mcManager->AddFile(tra3File);
+
+ run->AddTask(mcManager);
+ }
+ // ------------------------------------------------------------------------
+
+
+ // ----- Match reco to MC ------
+ if (bUseMC) {
+ CbmMatchRecoToMC* matchTask = new CbmMatchRecoToMC();
+ run->AddTask(matchTask);
+ }
+
+ // ----- Geometry interface initializer for tracking
+ run->AddTask(new CbmTrackingDetectorInterfaceInit());
+
+
+ // *******************************
+ // ** Detector QA tasks **
+ // *******************************
+
+ // ----- MUCH QA ---------------------------------
+ if (CbmSetup::Instance()->IsActive(ECbmModuleId::kMuch)) {
+ run->AddTask(new CbmMuchTransportQa());
+ run->AddTask(new CbmMuchDigitizerQa());
+ CbmMuchHitFinderQa* muchHitFinderQa = new CbmMuchHitFinderQa();
+ muchHitFinderQa->SetGeoFileName(muchParFile);
+ run->AddTask(muchHitFinderQa);
+ }
+
+ // ----- TRD QA ---------------------------------
+ if (CbmSetup::Instance()->IsActive(ECbmModuleId::kTrd)) {
+ run->AddTask(new CbmTrdMCQa());
+ //run->AddTask(new CbmTrdHitRateQa()); //opens lots of windows
+ //run->AddTask(new CbmTrdDigitizerPRFQa()); //works put currently doesn't do anything
+ //run->AddTask(new CbmTrdHitRateFastQa()); //opens lots of windows
+ run->AddTask(new CbmTrdHitProducerQa());
+ run->AddTask(new CbmTrdCalibTracker());
+ run->AddTask(new CbmTrackerInputQaTrd()); // Tracker requirements to TRD
+ }
+ // ------------------------------------------------------------------------
+
+ // ----- TRD QA ---------------------------------
+ if (CbmSetup::Instance()->IsActive(ECbmModuleId::kTof)) {
+ // TODO: SZh 19.10.2022:
+ // After the proper TOF digi parameters initialization it is appeared,
+ // that CbmTrackerInputQaTof causes segmentation violation in
+ // CbmTrackerInputQaTof::ResolutionQa() for the event-based running
+ // (before the branch with the related code inside the task was never
+ // executed). Temporarily commented that task, till the bug would not be
+ // resolved.
+
+ //run->AddTask(new CbmTrackerInputQaTof()); // Tracker requirements to TRD
+ }
+ // ------------------------------------------------------------------------
+
+ // ----- STS QA ---------------------------------
+ if (CbmSetup::Instance()->IsActive(ECbmModuleId::kSts)) {
+ //run->AddTask(new CbmStsDigitizeQa()); //opens lots of windows
+ run->AddTask(new CbmStsFindTracksQa());
+ run->AddTask(new CbmTrackingInputQaSts());
+ run->AddTask(new CbmCaInputQaSts(verbose, bUseMC));
+ }
+ // ------------------------------------------------------------------------
+
+ // ----- Event builder QA ---------------------------------
+ CbmBuildEventsQa* evBuildQA = new CbmBuildEventsQa();
+ run->AddTask(evBuildQA);
+ // ------------------------------------------------------------------------
+
+ // ----- Parameter database --------------------------------------------
+ std::cout << std::endl << std::endl;
+ std::cout << "-I- " << myName << ": Set runtime DB" << std::endl;
+ FairRuntimeDb* rtdb = run->GetRuntimeDb();
+ FairParRootFileIo* parIo1 = new FairParRootFileIo();
+ parIo1->open(parFile.Data(), "in");
+ rtdb->setFirstInput(parIo1);
+ if (!parFileList->IsEmpty()) {
+ FairParAsciiFileIo* parIo2 = new FairParAsciiFileIo();
+ parIo2->open(parFileList, "in");
+ rtdb->setSecondInput(parIo2);
+ }
+ rtdb->print();
+ // ------------------------------------------------------------------------
+
+ // ----- Run initialisation -------------------------------------------
+ std::cout << std::endl;
+ std::cout << "-I- " << myName << ": Initialise run" << std::endl;
+ run->Init();
+
+ // ----- Start run ----------------------------------------------------
+ std::cout << std::endl << std::endl;
+ std::cout << "-I- " << myName << ": Starting run" << std::endl;
+ run->Run(0, nEvents);
+ // ------------------------------------------------------------------------
+
+ // ----- Finish -------------------------------------------------------
+ timer.Stop();
+
+ FairMonitor::GetMonitor()->Print();
+ Double_t rtime = timer.RealTime();
+ Double_t ctime = timer.CpuTime();
+ std::cout << std::endl << std::endl;
+ std::cout << "Macro finished successfully." << std::endl;
+ std::cout << "Output file is " << sinkFile << std::endl;
+ std::cout << "Parameter file is " << parFile << std::endl;
+ std::cout << "Real time " << rtime << " s, CPU time " << ctime << " s" << std::endl;
+ std::cout << std::endl;
+ std::cout << " Test passed" << std::endl;
+ std::cout << " All ok " << std::endl;
+ // ------------------------------------------------------------------------
+
+ // ----- Resource monitoring ------------------------------------------
+ // Extract the maximal used memory an add is as Dart measurement
+ // This line is filtered by CTest and the value send to CDash
+ FairSystemInfo sysInfo;
+ Float_t maxMemory = sysInfo.GetMaxMemory();
+ std::cout << "<DartMeasurement name=\"MaxMemory\" type=\"numeric/double\">";
+ std::cout << maxMemory;
+ std::cout << "</DartMeasurement>" << std::endl;
+
+ Float_t cpuUsage = ctime / rtime;
+ std::cout << "<DartMeasurement name=\"CpuLoad\" type=\"numeric/double\">";
+ std::cout << cpuUsage;
+ std::cout << "</DartMeasurement>" << std::endl;
+ // ------------------------------------------------------------------------
+
+ // ----- Function needed for CTest runtime dependency -----------------
+ RemoveGeoManager();
+ // ------------------------------------------------------------------------
+}
diff --git a/reco/L1/CMakeLists.txt b/reco/L1/CMakeLists.txt
index df1bb9c9ab..e9c9c00351 100644
--- a/reco/L1/CMakeLists.txt
+++ b/reco/L1/CMakeLists.txt
@@ -77,6 +77,8 @@ set(SRCS
qa/CbmTrackerInputQaTrd.cxx
qa/CbmTrackerInputQaTof.cxx
qa/CbmTrackingInputQaSts.cxx
+ qa/CbmCaInputQaMuch.cxx
+ qa/CbmCaInputQaSts.cxx
)
set(NO_DICT_SRCS
@@ -95,6 +97,7 @@ set(HEADERS
L1Algo/L1Def.h
L1Algo/L1Vector.h
catools/CaToolsWindowFinder.h
+ qa/CbmCaQaCuts.h
)
diff --git a/reco/L1/CbmCaMCModule.h b/reco/L1/CbmCaMCModule.h
new file mode 100644
index 0000000000..1c7f4ea633
--- /dev/null
+++ b/reco/L1/CbmCaMCModule.h
@@ -0,0 +1,507 @@
+/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// \file CbmCaMCModule.h
+/// \brief CA Tracking performance interface for CBM (header)
+/// \since 23.09.2022
+/// \author S.Zharko <s.zharko@gsi.de>
+
+#ifndef CbmCaMCModule_h
+#define CbmCaMCModule_h 1
+
+#include "CbmL1DetectorID.h"
+#include "CbmL1Track.h"
+#include "CbmLink.h"
+#include "CbmMCDataArray.h"
+#include "CbmMCEventList.h"
+#include "CbmMCTrack.h"
+#include "CbmMatch.h"
+#include "CbmMuchPoint.h"
+#include "CbmMvdPoint.h"
+#include "CbmStsPoint.h"
+#include "CbmTimeSlice.h"
+#include "CbmTofPoint.h"
+#include "CbmTrdPoint.h"
+
+#include "TClonesArray.h"
+#include "TDatabasePDG.h"
+
+#include <numeric>
+#include <string>
+#include <string_view>
+#include <type_traits>
+
+#include "CaToolsMCData.h"
+#include "CaToolsMCPoint.h"
+#include "CaToolsQa.h"
+#include "L1Constants.h"
+#include "L1Parameters.h"
+#include "L1Undef.h"
+
+class CbmEvent;
+class CbmMCDataObject;
+class CbmL1Hit;
+class CbmL1Track;
+
+enum class L1DetectorID;
+
+/// Class CbmCaPerformcance is an interface to communicate between
+class CbmCaMCModule {
+public:
+ // *****************************************
+ // ** Constructors and destructor **
+ // *****************************************
+
+ /// Constructor
+ /// \param verbosity Verbosity level
+ CbmCaMCModule(int verb = 0, int perfMode = 3) : fVerbose(verb), fPerformanceMode(perfMode) {}
+
+ /// Destructor
+ ~CbmCaMCModule() = default;
+
+ /// Copy constructor
+ CbmCaMCModule(const CbmCaMCModule&) = delete;
+
+ /// Move constructor
+ CbmCaMCModule(CbmCaMCModule&&) = delete;
+
+ /// Copy assignment operator
+ CbmCaMCModule& operator=(const CbmCaMCModule&) = delete;
+
+ /// Move assignment operator
+ CbmCaMCModule& operator=(CbmCaMCModule&&) = delete;
+
+
+ // *****************************************
+ // ** Action definition functions **
+ // *****************************************
+
+ /// Defines performance action in the beginning of the run
+ /// \return Success flag
+ bool InitRun();
+
+ /// \brief Defines performance action in the beginning of each event or timeslice
+ /// \note This function should be called before hits initialization
+ /// \param pEvent Pointer to a current CbmEvent
+ void InitEvent(CbmEvent* pEvent);
+
+ /// \brief Processes event
+ /// Fills histograms and tables, should be called after the tracking done
+ void ProcessEvent(CbmEvent* pEvent);
+
+ /// \brief Initialize information about arrangement of points and hits of MC tracks within stations and the status
+ /// of track reconstructability
+ /// \param vHits Vector of hit objects
+ /// Calculates max number of points and hits on a station, number of consecutive stations containing a hit or point
+ /// and number of stations and points with hits. Provides the determination of track reconstructablity status
+ void InitTrackInfo(const L1Vector<CbmL1Hit>& vHits);
+
+ /// Defines performance action in the end of the run
+ void Finish();
+
+ /// Gets the first point index for a given detector subsystem
+ int GetFirstPointIndex(L1DetectorID detID) const
+ {
+ return std::accumulate(fvNofPointsUsed.cbegin(), fvNofPointsUsed.cbegin() + int(detID), 0);
+ }
+
+ /// Gets the first point index for a given detector subsystem
+ int GetLastPointIndex(L1DetectorID detID) const
+ {
+ return std::accumulate(fvNofPointsUsed.cbegin(), fvNofPointsUsed.cbegin() + int(detID) + 1, 0) - 1;
+ }
+
+ /// \brief Matches hit with MC point
+ /// \tparam DetId Detector ID
+ /// \param iHit External index of hit
+ /// \return MC-point index in fvMCPoints array
+ /// This method finds a match for a given hit or matches for hits clusters (in case of STS), finds the best
+ /// link in the match and returns the corresponding global ID of the MC points
+ template<L1DetectorID DetId>
+ int MatchHitWithMc(int iHit) const;
+
+ /// Assigns MC point indexes to each hit and hit indexes to each MC point
+ /// \param vInputExtHits Vector of external input hits
+ void MatchPointsWithHits(L1Vector<CbmL1Hit>& vInputExtHits);
+
+ /// \brief Matches reconstructed tracks with MC tracks
+ /// \param vRecoTracks Vector of reconstructed tracks
+ /// \param vInputExtHits Vector of external input hits
+ /// In the procedure, the maps of associated MC track indexes to corresponding number of hits are filled out and the
+ /// max purity is calculated for each reconstructed track in the TS/event. If the value of purity for a given MC track
+ /// is larger then a threshold, the corresponding track index is saved to the reconstructed track object, in the same
+ /// time the index of the reconstructed track object is saved to this MC track object. If purity for the MC track does
+ /// not pass the threshold, its index will not be accounted as a matched track, and the index of reconstructed track
+ /// will be added as an index of "touched" track.
+ void MatchRecoAndMCTracks(L1Vector<CbmL1Track>& vRecoTracks, L1Vector<CbmL1Hit>& vInputExtHits);
+
+
+ // ***********************
+ // ** Accessors **
+ // ***********************
+
+ /// Gets reference to MC data object
+ const ca::tools::MCData& GetMCData() const { return fMCData; }
+
+ /// Gets mutable pointer to MC data object
+ ca::tools::MCData* GetMCData() { return &fMCData; }
+
+ /// Registers pointer to the tracking parameters object
+ void SetParameters(const L1Parameters* pParameters) { fpParameters = pParameters; }
+
+ /// Sets used detector subsystems
+ /// \note Should be called before this->Init()
+ /// \param detID Id of detector
+ /// \param flag Flag: true - detector is used
+ void SetDetector(L1DetectorID detID, bool flag);
+
+ /// Sets legacy event mode:
+ /// \param flag: true - runs on events base, false - runs on time-slice base
+ void SetLegacyEventMode(bool flag) { fbLegacyEventMode = flag; }
+
+ /// Sets pointer to the vector of reconstructed tracks from reference
+ void SetRecoTrackContainer(const L1Vector<CbmL1Track>& vRecoTracks) { fpvRecoTracks = &vRecoTracks; }
+
+ /// Sets pointer to the vector of hits from reference
+ void SetHitContainer(const L1Vector<CbmL1Hit>& vHits) { fpvHits = &vHits; }
+
+private:
+ // ****************************************
+ // ** Local constant expressions **
+ // ****************************************
+
+ static constexpr int kNdetectors = L1Constants::size::kMaxNdetectors;
+
+ // *******************************
+ // ** Utility functions **
+ // *******************************
+
+ /// Calculates global index of MC point as a sum of a given local index and total provided number of points in previous
+ /// detector subsystem
+ /// \param iPointLocal Local index of MC point
+ /// \param detID Detector ID
+ int CalcGlobMCPointIndex(int iPointLocal, L1DetectorID detID) const
+ {
+ return iPointLocal + std::accumulate(fvNofPointsOrig.cbegin(), fvNofPointsOrig.cbegin() + int(detID), 0);
+ }
+
+ /// Checks class initialization. Throws std::logic_error, if initialization is incomplete at initialization call
+ void CheckInit() const;
+
+ /// Reads MC tracks from external trees and saves them to MCDataObject
+ void ReadMCTracks();
+
+ /// Reads MC points from external trees and saves them to MCDataObject
+ void ReadMCPoints();
+
+ /// Reads MC points in particular detector
+ template<L1DetectorID DetID>
+ void ReadMCPointsForDetector(CbmMCDataArray* pPoints);
+
+ /// Fills a single detector-specific MC point
+ /// \tparam DetID Detector subsystem ID
+ /// \param[in] iExtId Index of point in the external points container
+ /// \param[in] iEvent EventID of point in the external points container
+ /// \param[in] iFile FileID of point int the external points container
+ /// \param[out] intMCPoint Reference to the internal tracking MC point object (ca::tools::MCData)
+ /// \return Success flag:
+ /// - #true: Point is correct and is to be saved to the MCData object
+ /// - #false: Point is incorrect and will be ignored
+ template<L1DetectorID DetID>
+ bool FillMCPoint(int iExtId, int iEvent, int iFile, ca::tools::MCPoint& point);
+
+
+ // *******************************
+ // ** Utility variables **
+ // *******************************
+
+ const L1Parameters* fpParameters = nullptr; ///< Pointer to tracking parameters object
+
+ // ** Flags **
+
+ bool fbUseMvd = false; ///< is MVD used
+ bool fbUseSts = false; ///< is STS used
+ bool fbUseMuch = false; ///< is MuCh used
+ bool fbUseTrd = false; ///< is TRD used
+ bool fbUseTof = false; ///< is TOF used
+
+ bool fbLegacyEventMode = false; ///< if tracking uses events instead of time-slices (back compatibility)
+ int fVerbose = 0; ///< Verbosity level
+ int fPerformanceMode = -1; ///< Mode of performance
+
+ // *********************************
+ // ** Input data branches **
+ // *********************************
+
+ const CbmTimeSlice* fpTimeSlice = nullptr; ///< Current time slice
+
+ // ------ Mc-event
+ CbmMCEventList* fpMCEventList = nullptr; ///< MC event list
+ CbmMCDataObject* fpMCEventHeader = nullptr; ///< MC event header
+ CbmMCDataArray* fpMCTracks = nullptr; ///< MC tracks input
+
+ // ------ Mc-points
+ CbmMCDataArray* fpMvdPoints = nullptr; ///< MVD MC-points input container
+ CbmMCDataArray* fpStsPoints = nullptr; ///< STS MC-points input container
+ CbmMCDataArray* fpMuchPoints = nullptr; ///< MuCh MC-points input container
+ CbmMCDataArray* fpTrdPoints = nullptr; ///< TRD MC-points input container
+ CbmMCDataArray* fpTofPoints = nullptr; ///< TOF MC-points input container
+
+ std::array<int, kNdetectors> fvNofPointsOrig = {}; ///< Number of points by detector provided
+ std::array<int, kNdetectors> fvNofPointsUsed = {}; ///< Number of points used in performance
+
+ // ------ Matches
+ TClonesArray* fpMvdHitMatches = nullptr; ///< Array of MVD hit matches
+ TClonesArray* fpStsHitMatches = nullptr; ///< Array of STS hit matches
+ TClonesArray* fpMuchHitMatches = nullptr; ///< Array of MuCh hit matches
+ TClonesArray* fpTrdHitMatches = nullptr; ///< Array of TRD hit matches
+ TClonesArray* fpTofHitMatches = nullptr; ///< Array of TOF hit matches
+
+ TClonesArray* fpStsHits = nullptr; ///< Array of STS hits (currently needed for matching)
+ TClonesArray* fpStsClusterMatches = nullptr; ///< Array of STS cluster matches
+
+ // ------ Matching information
+ std::set<std::pair<int, int>> fFileEventIDs; ///< Set of file and event indexes: first - iFile, second - iEvent
+
+ // ************************************
+ // ** Current MC data object **
+ // ************************************
+
+ ca::tools::MCData fMCData {}; ///< MC-information in CA tracking internal format (tracks and points)
+
+ std::unique_ptr<ca::tools::Qa> fpQaModule = nullptr; ///< Pointer to CA QA module
+
+ // ********************************************
+ // ** Pointers to reconstructed data **
+ // ********************************************
+
+ const L1Vector<CbmL1Track>* fpvRecoTracks = nullptr; ///< Pointer to reconstructed track container
+ const L1Vector<CbmL1Hit>* fpvHits = nullptr; ///< Pointer to hits container
+};
+
+// **********************************************
+// ** Template function implementation **
+// **********************************************
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+template<L1DetectorID DetID>
+bool CbmCaMCModule::FillMCPoint(int iExtId, int iEvent, int iFile, ca::tools::MCPoint& point)
+{
+ // ----- Get positions, momenta, time and track ID
+ TVector3 posIn; // Position at entrance to station [cm]
+ TVector3 posOut; // Position at exist of station [cm]
+ TVector3 momIn; // 3-momentum at entrance to station [GeV/c]
+ TVector3 momOut; // 3-momentum at exit of station [GeV/c]
+ double time = undef::kD; // Time of MC point (after beginning of the event) [ns]
+ int iTmcExt = undef::kI32; // Track ID in the external container
+ // MVD
+ if constexpr (L1DetectorID::kMvd == DetID) {
+ auto* pExtPoint = L1_DYNAMIC_CAST<CbmMvdPoint*>(fpMvdPoints->Get(iFile, iEvent, iExtId));
+ if (!pExtPoint) {
+ LOG(warn) << "CbmCaMCModule: MVD MC point with iExtId = " << iExtId << ", iEvent = " << iEvent
+ << ", iFile = " << iFile << " does not exist";
+ return false;
+ }
+ pExtPoint->Position(posIn);
+ pExtPoint->PositionOut(posOut);
+ pExtPoint->Momentum(momIn);
+ pExtPoint->MomentumOut(momOut);
+ time = pExtPoint->GetTime();
+ iTmcExt = pExtPoint->GetTrackID();
+ }
+ // STS
+ else if constexpr (L1DetectorID::kSts == DetID) {
+ auto* pExtPoint = L1_DYNAMIC_CAST<CbmStsPoint*>(fpStsPoints->Get(iFile, iEvent, iExtId));
+ if (!pExtPoint) {
+ LOG(warn) << "CbmCaMCModule: STS MC point with iExtId = " << iExtId << ", iEvent = " << iEvent
+ << ", iFile = " << iFile << " does not exist";
+ return false;
+ }
+ pExtPoint->Position(posIn);
+ pExtPoint->PositionOut(posOut);
+ pExtPoint->Momentum(momIn);
+ pExtPoint->MomentumOut(momOut);
+ time = pExtPoint->GetTime();
+ iTmcExt = pExtPoint->GetTrackID();
+ }
+ // MuCh
+ else if constexpr (L1DetectorID::kMuch == DetID) {
+ auto* pExtPoint = L1_DYNAMIC_CAST<CbmMuchPoint*>(fpMuchPoints->Get(iFile, iEvent, iExtId));
+ if (!pExtPoint) {
+ LOG(warn) << "CbmCaMCModule: MuCh MC point with iExtId = " << iExtId << ", iEvent = " << iEvent
+ << ", iFile = " << iFile << " does not exist";
+ return false;
+ }
+ pExtPoint->Position(posIn);
+ pExtPoint->PositionOut(posOut);
+ pExtPoint->Momentum(momIn);
+ pExtPoint->Momentum(momOut);
+ time = pExtPoint->GetTime();
+ iTmcExt = pExtPoint->GetTrackID();
+ }
+ // TRD
+ else if constexpr (L1DetectorID::kTrd == DetID) {
+ auto* pExtPoint = L1_DYNAMIC_CAST<CbmTrdPoint*>(fpTrdPoints->Get(iFile, iEvent, iExtId));
+ if (!pExtPoint) {
+ LOG(warn) << "CbmCaMCModule: TRD MC point with iExtId = " << iExtId << ", iEvent = " << iEvent
+ << ", iFile = " << iFile << " does not exist";
+ return false;
+ }
+ pExtPoint->Position(posIn);
+ pExtPoint->PositionOut(posOut);
+ pExtPoint->Momentum(momIn);
+ pExtPoint->MomentumOut(momOut);
+ time = pExtPoint->GetTime();
+ iTmcExt = pExtPoint->GetTrackID();
+ }
+ // TOF
+ else if constexpr (L1DetectorID::kTof == DetID) {
+ auto* pExtPoint = L1_DYNAMIC_CAST<CbmTofPoint*>(fpTofPoints->Get(iFile, iEvent, iExtId));
+ if (!pExtPoint) {
+ LOG(warn) << "CbmCaMCModule: TOF MC point with iExtId = " << iExtId << ", iEvent = " << iEvent
+ << ", iFile = " << iFile << " does not exist";
+ return false;
+ }
+ pExtPoint->Position(posIn);
+ pExtPoint->Position(posOut);
+ pExtPoint->Momentum(momIn);
+ pExtPoint->Momentum(momOut);
+ time = pExtPoint->GetTime();
+ iTmcExt = pExtPoint->GetTrackID();
+ }
+ if (iTmcExt < 0) {
+ LOG(warn) << "CbmCaMCModule: For MC point with iExtId = " << iExtId << ", iEvent = " << iEvent
+ << ", iFile = " << iFile << " MC track is undefined (has ID = " << iTmcExt << ')';
+ return false;
+ }
+ TVector3 posMid = 0.5 * (posIn + posOut);
+ TVector3 momMid = 0.5 * (momIn + momOut);
+
+ // // ----- Get station index
+ int stationID = undef::kI32; // Index of station in the array of active tracking stations
+ int nStationsGeo = fpParameters->GetNstationsGeometry(DetID);
+ // MVD, STS TODO: Try to extend for MuCh and TRD
+ if constexpr (L1DetectorID::kMvd == DetID || L1DetectorID::kSts == DetID) {
+ // Determine active station global index
+ double bestDist = 1.e+20;
+ for (int iStLocGeo = 0; iStLocGeo < nStationsGeo; ++iStLocGeo) {
+ int iStActive = fpParameters->GetStationIndexActive(iStLocGeo, DetID);
+ if (iStActive < 0) { continue; }
+ const L1Station& station = fpParameters->GetStation(iStActive);
+ double dist = posIn.Z() - station.fZ[0];
+ if (fabs(dist) < fabs(bestDist)) {
+ bestDist = dist;
+ stationID = iStActive;
+ }
+ }
+ }
+ // MuCh, TRD
+ else if constexpr (L1DetectorID::kMuch == DetID || L1DetectorID::kTrd == DetID) {
+ // Offset for MuCh - 2.5 cm, for TRD - 20. cm. TODO: Where did these values came from?
+ double offset = (L1DetectorID::kMuch == DetID) ? 2.5 : 20.;
+ for (int iStLocGeo = 0; iStLocGeo < nStationsGeo; ++iStLocGeo) {
+ int iStActive = fpParameters->GetStationIndexActive(iStLocGeo, DetID);
+ if (iStActive < 0) { continue; }
+ const L1Station& station = fpParameters->GetStation(iStActive);
+ if (posMid.Z() > station.fZ[0] - offset) { stationID = iStActive; }
+ }
+ }
+
+ // ----- Reject MC points falling out of the time slice
+ // STS, MuCh, TRD, TOF
+ if constexpr (DetID != L1DetectorID::kMvd) {
+ if (!fbLegacyEventMode) { // fpTimeSlice != nullptr
+ double startT = fpTimeSlice->GetStartTime();
+ double endT = fpTimeSlice->GetEndTime();
+ double mcGlobT = time + fpMCEventList->GetEventTime(iEvent, iFile);
+
+ if ((startT > 0. && mcGlobT < startT) || (endT > 0. && mcGlobT > endT)) {
+ LOG(warn) << "CbmCaMCModule: MC point with iExtId = " << iExtId << ", iEvent = " << iEvent
+ << ", iFile = " << iFile << " and det id " << int(DetID) << " fell out of the TS duration [" << startT
+ << ", " << endT << "] with measured time = " << mcGlobT << " [ns]";
+ return false;
+ }
+ }
+ }
+
+ // ----- Fill MC point
+ point.SetEventId(iEvent);
+ point.SetFileId(iFile);
+ point.SetTime(time);
+
+ point.SetX(posMid.X());
+ point.SetY(posMid.Y());
+ point.SetZ(posMid.Z());
+
+ point.SetXIn(posIn.X());
+ point.SetYIn(posIn.Y());
+ point.SetZIn(posIn.Z());
+
+ point.SetXOut(posOut.X());
+ point.SetYOut(posOut.Y());
+ point.SetZOut(posOut.Z());
+
+ point.SetPx(momMid.X());
+ point.SetPy(momMid.Y());
+ point.SetPz(momMid.Z());
+
+ point.SetPxIn(momIn.X());
+ point.SetPyIn(momIn.Y());
+ point.SetPzIn(momIn.Z());
+
+ point.SetPxOut(momOut.X());
+ point.SetPyOut(momOut.Y());
+ point.SetPzOut(momOut.Z());
+
+ int iTmcInt = fMCData.FindInternalTrackIndex(iTmcExt, iEvent, iFile);
+
+ point.SetId(fMCData.GetNofPoints()); // select current number of points as a local id of points
+ point.SetTrackId(iTmcInt);
+ point.SetStationId(stationID);
+ point.SetDetectorId(DetID);
+
+ auto* pExtTrk = L1_DYNAMIC_CAST<CbmMCTrack*>(fpMCTracks->Get(iFile, iEvent, iTmcExt));
+ if (!pExtTrk) {
+ LOG(warn) << "CbmCaMCModule: MC track with iTmcExt = " << iTmcExt << ", iEvent = " << iEvent
+ << ", iFile = " << iFile << " MC track is undefined (nullptr)";
+ return false;
+ }
+ point.SetPdgCode(pExtTrk->GetPdgCode());
+ point.SetMotherId(pExtTrk->GetMotherId());
+
+ auto* pPdgDB = TDatabasePDG::Instance()->GetParticle(point.GetPdgCode());
+ point.SetMass(pPdgDB ? pPdgDB->Mass() : 0.); /// TODO: Set from track
+ point.SetCharge(pPdgDB ? pPdgDB->Charge() / 3. : 0.);
+
+ return true;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+// NOTE: template is used, because another template function FillMCPoint is used inside
+template<L1DetectorID DetID>
+void CbmCaMCModule::ReadMCPointsForDetector(CbmMCDataArray* pPoints)
+{
+ if constexpr (L1DetectorID::kTof != DetID) {
+ fvNofPointsUsed[int(DetID)] = 0;
+ for (const auto& key : fFileEventIDs) {
+ int iFile = key.first;
+ int iEvent = key.second;
+ int nPointsEvent = pPoints->Size(iFile, iEvent);
+ for (int iP = 0; iP < nPointsEvent; ++iP) {
+ ca::tools::MCPoint aPoint;
+ if (FillMCPoint<DetID>(iP, iEvent, iFile, aPoint)) {
+ aPoint.SetExternalId(CalcGlobMCPointIndex(iP, DetID));
+ int iPInt = fMCData.GetNofPoints(); // assign an index of point in internal container
+ if (aPoint.GetTrackId() > -1) { fMCData.GetTrack(aPoint.GetTrackId()).AddPointIndex(iPInt); }
+ fMCData.AddPoint(aPoint);
+ ++fvNofPointsUsed[int(DetID)];
+ }
+ } // iP: end
+ } // key: end
+ }
+}
+
+
+#endif // CbmCaMCModule_h
diff --git a/reco/L1/L1LinkDef.h b/reco/L1/L1LinkDef.h
index 7a1cccc9b4..ad266732d7 100644
--- a/reco/L1/L1LinkDef.h
+++ b/reco/L1/L1LinkDef.h
@@ -30,6 +30,8 @@
#pragma link C++ class CbmTrackingInputQaSts + ;
#pragma link C++ class CbmTrackerInputQaTrd + ;
#pragma link C++ class CbmTrackerInputQaTof + ;
+#pragma link C++ class CbmCaInputQaMuch + ;
+#pragma link C++ class CbmCaInputQaSts + ;
#pragma link C++ class ca::tools::WindowFinder + ;
#endif
diff --git a/reco/L1/qa/CbmCaInputQaSts.cxx b/reco/L1/qa/CbmCaInputQaSts.cxx
new file mode 100644
index 0000000000..7193fd7f17
--- /dev/null
+++ b/reco/L1/qa/CbmCaInputQaSts.cxx
@@ -0,0 +1,1193 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// \file CbmCaInputQaSts.cxx
+/// \date 13.01.2023
+/// \brief QA-task for CA tracking input from MuCh detector (implementation)
+/// \author S.Zharko <s.zharko@gsi.de>
+
+#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 "CbmStsCluster.h"
+#include "CbmStsHit.h"
+#include "CbmStsPoint.h"
+#include "CbmStsTrackingInterface.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(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()
+{
+ 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
+{
+ LOG_IF(info, fVerbose > 0) << fName << ": Checking QA ...";
+
+ bool res = true;
+
+ int nSt = fpDetInterface->GetNtrackingStations();
+
+
+ // **************************************************************
+ // ** Basic checks, available both for real and simulated data **
+ // **************************************************************
+
+ // ----- Checks for mismatches in the order of 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(-cuts::kMaxDzStHitSts);
+ int iBinMax = fvph_hit_station_delta_z[iSt]->FindBin(+cuts::kMaxDzStHitSts);
+
+ 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
+ //
+ 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.";
+ 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;
+ }
+ }
+
+ // ----- 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
+
+ // 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();
+ };
+
+ 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)";
+ 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;
+ }
+ }
+
+ LOG_IF(error, !res) << fName << ": QA check failed";
+
+ return res;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void CbmCaInputQaSts::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_du.clear();
+ fvph_hit_dv.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_u.clear();
+ fvph_res_v.clear();
+ fvph_res_t.clear();
+
+ fvph_pull_x.clear();
+ fvph_pull_y.clear();
+ fvph_pull_u.clear();
+ fvph_pull_v.clear();
+ fvph_pull_t.clear();
+
+ fvph_res_x_vs_x.clear();
+ fvph_res_y_vs_y.clear();
+ fvph_res_u_vs_u.clear();
+ fvph_res_v_vs_v.clear();
+ fvph_res_t_vs_t.clear();
+
+ fvph_pull_x_vs_x.clear();
+ fvph_pull_y_vs_y.clear();
+ fvph_pull_u_vs_u.clear();
+ fvph_pull_v_vs_v.clear();
+ fvph_pull_t_vs_t.clear();
+
+ fvpe_reco_eff_vs_xy.clear();
+ fvpe_reco_eff_vs_r.clear();
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void CbmCaInputQaSts::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 CbmStsHit*>(fpHits->At(iH));
+ LOG_IF(fatal, !pHit) << fName << ": hit with iH = " << iH << " is not an CbmStsHit (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;
+ double stPhiF = fpDetInterface->GetStripsStereoAngleFront(iSt); // STS front strips stereo angle
+ double stPhiB = fpDetInterface->GetStripsStereoAngleBack(iSt); // STS back strips stereo angle
+
+ // Hit position
+ double xHit = pHit->GetX();
+ double yHit = pHit->GetY();
+ double zHit = pHit->GetZ();
+ double uHit = xHit * cos(stPhiF) + yHit * sin(stPhiB);
+ double vHit = xHit * cos(stPhiB) + yHit * sin(stPhiB);
+
+ double duHit = pHit->GetDu();
+ double dvHit = pHit->GetDv();
+ 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_du[iSt]->Fill(duHit);
+ fvph_hit_dv[iSt]->Fill(dvHit);
+ 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_du[nSt]->Fill(duHit);
+ fvph_hit_dv[nSt]->Fill(dvHit);
+ fvph_hit_dt[nSt]->Fill(dtHit);
+
+
+ // **********************
+ // ** MC information QA
+
+ if (IsMCUsed()) {
+ CbmMatch hitMatch = GetHitMatch(iH);
+
+ // Evaluate number of hits per one MC point
+ int nMCpoints = 0; // Number of MC points for this hit
+ for (int iLink = 0; iLink < hitMatch.GetNofLinks(); ++iLink) {
+ const auto& link = hitMatch.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 = hitMatch.GetMatchedLink();
+
+ // Skip noisy links
+ if (bestPointLink.GetIndex() < 0) { continue; }
+
+ // Point matched by the best link
+ const auto* pMCPoint = dynamic_cast<const CbmStsPoint*>(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 uMC = xMC * cos(stPhiF) + yMC * sin(stPhiF);
+ //double vMC = xMC * cos(stPhiB) + yMC * sin(stPhiB);
+ 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 uMCs = xMCs * cos(stPhiF) + yMCs * sin(stPhiF);
+ double vMCs = xMCs * cos(stPhiB) + yMCs * sin(stPhiB);
+ double tMCs = tMC + shiftZ / (pzMC * phys::kSpeedOfLight) * sqrt(mass * mass + pMC * pMC);
+
+ // Residuals
+ double xRes = xHit - xMCs;
+ double yRes = yHit - yMCs;
+ double uRes = uHit - uMCs;
+ double vRes = vHit - vMCs;
+ double tRes = tHit - tMCs;
+
+ // ------ Cuts
+
+ if (std::fabs(pMCPoint->GetPzOut()) < 0.01) { 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_u[iSt]->Fill(uRes);
+ fvph_res_v[iSt]->Fill(vRes);
+ fvph_res_t[iSt]->Fill(tRes);
+
+ fvph_pull_x[iSt]->Fill(xRes / dxHit);
+ fvph_pull_y[iSt]->Fill(yRes / dyHit);
+ fvph_pull_u[iSt]->Fill(uRes / duHit);
+ fvph_pull_v[iSt]->Fill(vRes / dvHit);
+ 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_u_vs_u[iSt]->Fill(uHit, uRes);
+ fvph_res_v_vs_v[iSt]->Fill(vHit, vRes);
+ 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_u_vs_u[iSt]->Fill(uHit, uRes / duHit);
+ fvph_pull_v_vs_v[iSt]->Fill(vHit, vRes / dvHit);
+ 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_u[nSt]->Fill(uRes);
+ fvph_res_v[nSt]->Fill(vRes);
+ fvph_res_t[nSt]->Fill(tRes);
+
+ fvph_pull_x[nSt]->Fill(xRes / dxHit);
+ fvph_pull_y[nSt]->Fill(yRes / dyHit);
+ fvph_pull_u[nSt]->Fill(uRes / duHit);
+ fvph_pull_v[nSt]->Fill(vRes / dvHit);
+ 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_u_vs_u[nSt]->Fill(uHit, uRes);
+ fvph_res_v_vs_v[nSt]->Fill(vHit, vRes);
+ 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_u_vs_u[nSt]->Fill(uHit, uRes / duHit);
+ fvph_pull_v_vs_v[nSt]->Fill(vHit, vRes / dvHit);
+ 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 CbmStsPoint*>(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 isOneHitForOnePoint = (vNofHitsPerPoint[iE][iP] == 1);
+ bool isSevHitsForOnePoint = (vNofHitsPerPoint[iE][iP] > 1);
+
+ fvpe_reco_eff_vs_xy[iSt]->Fill(isOneHitForOnePoint, xMC, yMC);
+ fvpe_reco_eff_vs_xy[nSt]->Fill(isOneHitForOnePoint, xMC, yMC);
+
+ fvpe_reco_eff_vs_r[iSt]->Fill(isOneHitForOnePoint, rMC);
+ fvpe_reco_eff_vs_r[nSt]->Fill(isOneHitForOnePoint, rMC);
+
+ } // loop over MC-points: end
+ } // loop over MC-events: end
+ } // MC is used: end
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+InitStatus CbmCaInputQaSts::InitDataBranches()
+{
+ // STS tracking detector interface
+ fpDetInterface = CbmStsTrackingInterface::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("StsHit"));
+ LOG_IF(fatal, !fpHits) << "\033[1;31m" << fName << ": container of reconstructed hits in STS is not found\033[0m";
+
+ // Clusters container
+ fpClusters = dynamic_cast<TClonesArray*>(pFairRootManager->GetObject("StsCluster"));
+ LOG_IF(fatal, !fpClusters) << "\033[1;31m" << fName << ": container of hit clusters in STS 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("StsPoint");
+ LOG_IF(fatal, !fpMCTracks) << "\033[1;31m" << fName << ": MC point branch is not found for STS\033[0m";
+
+ // Cluster matches
+ fpClusterMatches = dynamic_cast<TClonesArray*>(pFairRootManager->GetObject("StsClusterMatch"));
+ LOG_IF(fatal, !fpClusterMatches) << "\033[1;31m]" << fName << ": cluster match branch is not found for STS\033[0m";
+ }
+
+ return kSUCCESS;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+InitStatus CbmCaInputQaSts::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("", "");
+ }
+
+ // x-coordinate
+ auto* pc_res_u = MakeCanvas<CbmQaCanvas>("res_u", "Residuals for u coordinate", 1600, 800);
+ pc_res_u->Divide2D(nSt);
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ pc_res_u->cd(iSt + 1);
+ fvph_res_u[iSt]->DrawCopy("", "");
+ }
+
+ // x-coordinate
+ auto* pc_res_v = MakeCanvas<CbmQaCanvas>("res_v", "Residuals for v coordinate", 1600, 800);
+ pc_res_v->Divide2D(nSt);
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ pc_res_v->cd(iSt + 1);
+ fvph_res_v[iSt]->DrawCopy("", "");
+ }
+
+ // x-coordinate
+ 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("", "");
+ }
+
+ // x-coordinate
+ auto* pc_pull_u = MakeCanvas<CbmQaCanvas>("pull_u", "Pulls for u coordinate", 1600, 800);
+ pc_pull_u->Divide2D(nSt);
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ pc_pull_u->cd(iSt + 1);
+ fvph_pull_u[iSt]->DrawCopy("", "");
+ }
+
+ // x-coordinate
+ auto* pc_pull_v = MakeCanvas<CbmQaCanvas>("pull_v", "Pulls for v coordinate", 1600, 800);
+ pc_pull_v->Divide2D(nSt);
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ pc_pull_v->cd(iSt + 1);
+ fvph_pull_v[iSt]->DrawCopy("", "");
+ }
+
+ // x-coordinate
+ 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 CbmCaInputQaSts::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);
+ fvph_hit_dv.resize(nSt + 1, nullptr);
+ fvph_hit_du.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 STS 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_du" + nsuff;
+ sT = (TString) "Hit position error across front strips" + tsuff + ";du_{hit} [cm]";
+ fvph_hit_du[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, kRHitDu[0], kRHitDu[1]);
+
+ sN = (TString) "hit_dv" + nsuff;
+ sT = (TString) "Hit position error across back strips" + tsuff + ";dv_{hit} [cm]";
+ fvph_hit_dv[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, kRHitDv[0], kRHitDv[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_u.resize(nSt + 1, nullptr);
+ fvph_res_v.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_u.resize(nSt + 1, nullptr);
+ fvph_pull_v.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_u_vs_u.resize(nSt + 1, nullptr);
+ fvph_res_v_vs_v.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_u_vs_u.resize(nSt + 1, nullptr);
+ fvph_pull_v_vs_v.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 STS 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 STS 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_u" + nsuff;
+ sT = (TString) "Residuals for front strips coordinate" + tsuff + ";u_{reco} - u_{MC} [cm]";
+ fvph_res_u[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, kRResU[0], kRResU[1]);
+
+ sN = (TString) "res_v" + nsuff;
+ sT = (TString) "Residuals for back strips coordinate" + tsuff + ";v_{reco} - v_{MC} [cm]";
+ fvph_res_v[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, kRResV[1], kRResV[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_u" + nsuff;
+ sT = (TString) "Pulls for front strips coordinate" + tsuff + ";(u_{reco} - u_{MC}) / #sigma_{u}^{reco}";
+ fvph_pull_u[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, kRPullU[0], kRPullU[1]);
+
+ sN = (TString) "pull_v" + nsuff;
+ sT = (TString) "Pulls for back strips coordinate" + tsuff + ";(v_{reco} - v_{MC}) / #sigma_{v}^{reco}";
+ fvph_pull_v[iSt] = MakeHisto<TH1F>(sN, sT, kNbins, kRPullV[0], kRPullV[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_u_vs_u" + nsuff;
+ sT = (TString) "Residuals for front strips coordinate" + tsuff + ";u_{reco} [cm];u_{reco} - u_{MC} [cm]";
+ fvph_res_u_vs_u[iSt] = MakeHisto<TH2F>(sN, sT, kNbins, 0, 0, kNbins, kRResU[0], kRResU[1]);
+
+ sN = (TString) "res_v_vs_u" + nsuff;
+ sT = (TString) "Residuals for back strips coordinate" + tsuff + ";v_{reco} [cm];v_{reco} - v_{MC} [cm]";
+ fvph_res_v_vs_v[iSt] = MakeHisto<TH2F>(sN, sT, kNbins, 0, 0, kNbins, kRResV[0], kRResV[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_u_vs_u" + nsuff;
+ sT = (TString) "Pulls for front strips coordinate" + tsuff
+ + ";u_{reco} [cm];(u_{reco} - u_{MC}) / #sigma_{u}^{reco}";
+ fvph_pull_u_vs_u[iSt] = MakeHisto<TH2F>(sN, sT, kNbins, 0, 0, kNbins, kRPullU[0], kRPullU[1]);
+
+ sN = (TString) "pull_v_vs_v" + nsuff;
+ sT =
+ (TString) "Pulls for back strips coordinate" + tsuff + ";v_{reco} [cm];(v_{reco} - v_{MC}) / #sigma_{v}^{reco}";
+ fvph_pull_v_vs_v[iSt] = MakeHisto<TH2F>(sN, sT, kNbins, 0, 0, kNbins, kRPullV[0], kRPullV[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<TEfficiency>(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);
+ }
+ }
+ return kSUCCESS;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+CbmMatch CbmCaInputQaSts::GetHitMatch(int iHit) const
+{
+ CbmMatch aMatch;
+
+ const auto* pHit = static_cast<const CbmStsHit*>(fpHits->At(iHit)); // NOTE: dynamic cast check is done in main loop
+
+ // Get front and back clusters
+ int iCf = pHit->GetFrontClusterId();
+ int iCb = pHit->GetBackClusterId();
+
+ // TODO: errors instead of fatals?
+ LOG_IF(fatal, iCf < 0) << fName << ": hit with id = " << iHit << " has wrong front cluster (iCf = " << iCf << ')';
+ LOG_IF(fatal, iCb < 0) << fName << ": hit with id = " << iHit << " has wrong back cluster (iCb = " << iCb << ')';
+
+ const auto* pClusterF = dynamic_cast<const CbmStsCluster*>(fpClusters->At(iCf));
+ const auto* pClusterB = dynamic_cast<const CbmStsCluster*>(fpClusters->At(iCb));
+
+ LOG_IF(fatal, !pClusterF) << fName << ": front cluster does not exist for hit with id = " << iHit
+ << " and cluster id = " << iCf;
+ LOG_IF(fatal, !pClusterB) << fName << ": back cluster does not exist for hit with id = " << iHit
+ << " and cluster id = " << iCb;
+
+ const auto* pClusterMatchF = dynamic_cast<const CbmMatch*>(fpClusterMatches->At(iCf));
+ const auto* pClusterMatchB = dynamic_cast<const CbmMatch*>(fpClusterMatches->At(iCb));
+
+ LOG_IF(fatal, !pClusterMatchF) << fName << ": front cluster match was not found for cluster with id = " << iCf;
+ LOG_IF(fatal, !pClusterMatchB) << fName << ": back cluster match was not found for cluster with id = " << iCf;
+
+ // Check addresses of hit and cluster
+ int addrHit = pHit->GetAddress();
+ int addrClusterF = pClusterF->GetAddress();
+ int addrClusterB = pClusterB->GetAddress();
+ bool ifAddrCons = addrHit == addrClusterF && addrClusterF == addrClusterB;
+ LOG_IF(fatal, !ifAddrCons) << fName << ": hit, front and (or) back cluster has inconsistent addresses: hit - "
+ << addrHit << ", front - " << addrClusterF << ", back - " << addrClusterB << '\n';
+
+ // CUSTOM MATCHING IN STS: The link is selected if it is presented BOTH in front and back cluster, thus only true
+ // hits are matched. The standard matching in STS assumes only one (front or back) cluster
+ // to be matched, so it matches both true and fake hits.
+
+ for (int iLinkF = 0; iLinkF < pClusterMatchF->GetNofLinks(); ++iLinkF) {
+ const auto& linkF = pClusterMatchF->GetLink(iLinkF);
+ for (int iLinkB = 0; iLinkB < pClusterMatchB->GetNofLinks(); ++iLinkB) {
+ const auto& linkB = pClusterMatchB->GetLink(iLinkB);
+ if (linkF == linkB) {
+ aMatch.AddLink(linkF);
+ aMatch.AddLink(linkB);
+ }
+ }
+ }
+ return aMatch;
+}
diff --git a/reco/L1/qa/CbmCaInputQaSts.h b/reco/L1/qa/CbmCaInputQaSts.h
new file mode 100644
index 0000000000..3e205c7803
--- /dev/null
+++ b/reco/L1/qa/CbmCaInputQaSts.h
@@ -0,0 +1,201 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// \file CbmCaInputQaSts.h
+/// \date 13.01.2023
+/// \brief QA-task for CA tracking input from MuCh detector (header)
+/// \author S.Zharko <s.zharko@gsi.de>
+
+
+#ifndef CbmCaInputQaSts_h
+#define CbmCaInputQaSts_h 1
+
+#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 CbmStsHit;
+class CbmStsTrackingInterface;
+class TClonesArray;
+class TH1F;
+class TH2F;
+class TEfficiency;
+
+/// A QA-task class, which provides assurance of MuCh hits and geometry
+class CbmCaInputQaSts : public CbmQaTask {
+public:
+ /// Constructor from parameters
+ /// \param verbose Verbose level
+ /// \param isMCUsed Flag, whether MC information is available for this task
+ CbmCaInputQaSts(int verbose, bool isMCUsed);
+
+ ClassDef(CbmCaInputQaSts, 0);
+
+protected:
+ // ********************************************
+ // ** Virtual method override from CbmQaTask **
+ // ********************************************
+
+ /// Checks results of the QA and returns some flag
+ bool Check() const;
+
+ /// Initializes data branches
+ InitStatus InitDataBranches();
+
+ /// Initializes canvases
+ InitStatus InitCanvases();
+
+ /// Initializes histograms
+ InitStatus InitHistograms();
+
+ /// 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 **
+ // *********************
+
+ /// Gets match object for the hit
+ /// \param iHit Index of hit
+ /// \return Match object
+ CbmMatch GetHitMatch(int iHit) const;
+
+
+ // *********************
+ // ** Class variables **
+ // *********************
+
+
+ // ** Data branches **
+ CbmStsTrackingInterface* fpDetInterface = nullptr; ///< Instance of detector interface
+
+ CbmTimeSlice* fpTimeSlice = nullptr; ///< Pointer to current time-slice
+
+ TClonesArray* fpHits = nullptr; ///< Array of hits
+ TClonesArray* fpClusters = nullptr; ///< Array of hit clusters
+
+ 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* 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]
+
+ 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]
+
+ 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] = {-.005, .005}; ///< Range for hit x coordinate [cm]
+ static constexpr double kRHitDy[2] = {-.005, .005}; ///< 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]
+
+ 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 kRResU[2] = {-5., 5.}; ///< Range for residual of v coordinate [cm]
+ static constexpr double kRResV[2] = {-.4, .4}; ///< 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 kRPullU[2] = {-400., 400.}; ///< Range for pull of v coordinate
+ 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
+ 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_du;
+ std::vector<TH1F*> fvph_hit_dv;
+ 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_u; ///< residual for u coordinate in different stations
+ std::vector<TH1F*> fvph_res_v; ///< residual for v 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_u_vs_u; ///< residual for u coordinate in different stations
+ std::vector<TH2F*> fvph_res_v_vs_v; ///< residual for v 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_u; ///< pull for u coordinate in different stations
+ std::vector<TH1F*> fvph_pull_v; ///< pull for v 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_u_vs_u; ///< pull for u coordinate in different stations
+ std::vector<TH2F*> fvph_pull_v_vs_v; ///< pull for v coordinate in different stations
+ 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);
+ }
+};
+
+#endif // CbmCaInputQaSts_h
diff --git a/reco/L1/qa/CbmCaQaCuts.h b/reco/L1/qa/CbmCaQaCuts.h
new file mode 100644
index 0000000000..08810ac00a
--- /dev/null
+++ b/reco/L1/qa/CbmCaQaCuts.h
@@ -0,0 +1,17 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+#ifndef CbmCaQaCuts_h
+#define CbmCaQaCuts_h 1
+
+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)
+
+ // Max difference between hit and station z
+ constexpr double kMaxDzStHitSts = 1.; ///< Max distance for STS
+}; // namespace cbm::ca::qa::cuts
+
+#endif // CbmCaQaCuts_h
diff --git a/reco/L1/qa/CbmTrackerInputQaTrd.cxx b/reco/L1/qa/CbmTrackerInputQaTrd.cxx
index 5ebe1fae51..04f61b5112 100644
--- a/reco/L1/qa/CbmTrackerInputQaTrd.cxx
+++ b/reco/L1/qa/CbmTrackerInputQaTrd.cxx
@@ -123,6 +123,8 @@ void CbmTrackerInputQaTrd::DeInit()
void CbmTrackerInputQaTrd::SetParContainers()
{
+ // FIXME: SZh 13.01.2023: With tracking detector interfaces these data bases are no longer needed
+ // to be initialized in this class
fTrdDigiPar = nullptr;
fTrdGeoPar = nullptr;
diff --git a/reco/L1/qa/CbmTrackingInputQaSts.cxx b/reco/L1/qa/CbmTrackingInputQaSts.cxx
index b1d32e068c..b05d2dcbcf 100644
--- a/reco/L1/qa/CbmTrackingInputQaSts.cxx
+++ b/reco/L1/qa/CbmTrackingInputQaSts.cxx
@@ -542,26 +542,11 @@ CbmMatch CbmTrackingInputQaSts::MatchHits(const CbmStsHit* pHit, int iHit)
//
double CbmTrackingInputQaSts::ParticleMass(int pdg)
{
+ // FIXME: SZh: Use masses defined by CbmMCTrack class instead
if (fabs(pdg) < 9999999 && ((TParticlePDG*) TDatabasePDG::Instance()->GetParticle(pdg))) {
return TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
}
- constexpr double kAmuToGev {0.93149410242}; // 1 amu in GeV/c2
-
- // Define masses of several light nuclei by hands
- if (fabs(pdg) == 1000010020) {
- return kAmuToGev * 2.01410177812; // deutron
- }
- else if (fabs(pdg) == 1000010030) {
- return kAmuToGev * 3.01604927790; // triton
- }
- else if (fabs(pdg) == 1000020030) {
- return kAmuToGev * 3.01602932007; // He-3
- }
- else if (fabs(pdg) == 1000020040) {
- return kAmuToGev * 4.00260325413; // He-4
- }
-
LOG(warn) << "\033[1;31m Found mass for pdg = " << pdg << " is undefined. "
<< "Please, provide data for this particle\033[0m\n";
@@ -753,8 +738,9 @@ void CbmTrackingInputQaSts::ResolutionQa()
double mass = ParticleMass(pdgCode);
constexpr double speedOfLight {29.9792458}; // cm/ns
TVector3 mom;
+ double mom2 = mcPx * mcPx + mcPy * mcPy + mcPz * mcPz;
- mcT += dz / (mcPz * speedOfLight) * sqrt(mass * mass + mom.Mag2());
+ mcT += dz / (mcPz * speedOfLight) * sqrt(mass * mass + mom2);
double dx = pHit->GetX() - mcX;
double dy = pHit->GetY() - mcY;
--
GitLab