From c8ace4dc9c65c5e9fa129c5b77c0745c8bf77c23 Mon Sep 17 00:00:00 2001
From: "s.zharko@gsi.de" <s.zharko@gsi.de>
Date: Fri, 15 Jul 2022 21:03:38 +0200
Subject: [PATCH] L1: Updates for STS tracking input QA
---
core/base/CbmTrackingDetectorInterfaceBase.h | 5 +
.../detectors/much/CbmMuchTrackingInterface.h | 12 +-
core/detectors/sts/CbmStsTrackingInterface.h | 10 +-
core/detectors/tof/CbmTofTrackingInterface.h | 13 +
core/detectors/trd/CbmTrdTrackingInterface.h | 10 +-
core/qa/CbmQaHist.h | 2 +-
mvd/CbmMvdTrackingInterface.h | 11 +-
reco/L1/CMakeLists.txt | 2 -
reco/L1/CbmL1ReadEvent.cxx | 4 +-
reco/L1/L1LinkDef.h | 1 -
reco/L1/qa/CbmTrackingInputQaBase.cxx | 319 -------
reco/L1/qa/CbmTrackingInputQaBase.h | 195 -----
reco/L1/qa/CbmTrackingInputQaSts.cxx | 787 +++++++++++++++++-
reco/L1/qa/CbmTrackingInputQaSts.h | 234 +++++-
14 files changed, 1031 insertions(+), 574 deletions(-)
delete mode 100644 reco/L1/qa/CbmTrackingInputQaBase.cxx
delete mode 100644 reco/L1/qa/CbmTrackingInputQaBase.h
diff --git a/core/base/CbmTrackingDetectorInterfaceBase.h b/core/base/CbmTrackingDetectorInterfaceBase.h
index dc9580829e..436c0a88b8 100644
--- a/core/base/CbmTrackingDetectorInterfaceBase.h
+++ b/core/base/CbmTrackingDetectorInterfaceBase.h
@@ -74,6 +74,11 @@ public:
/// \return Local index of the tracking station
virtual int GetTrackingStationIndex(const CbmPixelHit* hit) const = 0;
+ /// Gets a tracking station by the address of element
+ /// \param address Unique element address
+ /// \return Local index of the tracking station
+ virtual int GetTrackingStationIndex(int address) const = 0;
+
/// Gets max size of a station along the X-axis
/// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
/// \return Size of station along the X-axis [cm]
diff --git a/core/detectors/much/CbmMuchTrackingInterface.h b/core/detectors/much/CbmMuchTrackingInterface.h
index 0b637738a6..362288aeef 100644
--- a/core/detectors/much/CbmMuchTrackingInterface.h
+++ b/core/detectors/much/CbmMuchTrackingInterface.h
@@ -104,8 +104,16 @@ public:
/// \return Local index of the tracking station
int GetTrackingStationIndex(const CbmPixelHit* hit) const
{
- return CbmMuchGeoScheme::GetStationIndex(hit->GetAddress()) * 3
- + CbmMuchGeoScheme::GetLayerIndex(hit->GetAddress());
+ return GetTrackingStationIndex(hit->GetAddress());
+ }
+
+ /// Gets a tracking station by the address of element
+ /// \param address Unique element address
+ /// \return Local index of the tracking station
+ int GetTrackingStationIndex(int address) const
+ {
+ return CbmMuchGeoScheme::GetStationIndex(address) * 3
+ + CbmMuchGeoScheme::GetLayerIndex(address);
}
/// Gets max size of a station along the X-axis
diff --git a/core/detectors/sts/CbmStsTrackingInterface.h b/core/detectors/sts/CbmStsTrackingInterface.h
index d3ffb73619..d279ba6d73 100644
--- a/core/detectors/sts/CbmStsTrackingInterface.h
+++ b/core/detectors/sts/CbmStsTrackingInterface.h
@@ -109,7 +109,15 @@ public:
/// \return Local index of the tracking station
int GetTrackingStationIndex(const CbmPixelHit* hit) const
{
- return CbmStsSetup::Instance()->GetStationNumber(hit->GetAddress());
+ return GetTrackingStationIndex(hit->GetAddress());
+ }
+
+ /// Gets a tracking station by the address of element
+ /// \param address Unique element address
+ /// \return Local index of the tracking station
+ int GetTrackingStationIndex(int address) const
+ {
+ return CbmStsSetup::Instance()->GetStationNumber(address);
}
/// Gets max size of a station along the X-axis
diff --git a/core/detectors/tof/CbmTofTrackingInterface.h b/core/detectors/tof/CbmTofTrackingInterface.h
index 2c3f2e222a..1a0fa5fdba 100644
--- a/core/detectors/tof/CbmTofTrackingInterface.h
+++ b/core/detectors/tof/CbmTofTrackingInterface.h
@@ -111,6 +111,19 @@ public:
return iSt;
}
+ /// Gets a tracking station by the address of element
+ /// \param address Unique element address
+ /// \return Local index of the tracking station
+ int GetTrackingStationIndex(int /*address*/) const
+ {
+ // int iSt = fDigiBdfPar->GetTrackingStation(CbmTofAddress::GetSmType(address),
+ // CbmTofAddress::GetSmId(address),
+ // CbmTofAddress::GetRpcId(address));
+ // NOTE: Implement, when the "mcbm_beam_2021_07_surveyed" parameters will be fixed
+
+ return -1; // iSt;
+ }
+
/// Gets max size of a station along the X-axis
/// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
diff --git a/core/detectors/trd/CbmTrdTrackingInterface.h b/core/detectors/trd/CbmTrdTrackingInterface.h
index f193db4338..21108d9d0b 100644
--- a/core/detectors/trd/CbmTrdTrackingInterface.h
+++ b/core/detectors/trd/CbmTrdTrackingInterface.h
@@ -98,7 +98,15 @@ public:
/// Gets a tracking station of a CbmPixelHit
/// \param hit A pointer to CbmPixelHit
/// \return Local index of the tracking station
- int GetTrackingStationIndex(const CbmPixelHit* hit) const { return CbmTrdAddress::GetLayerId(hit->GetAddress()); }
+ int GetTrackingStationIndex(const CbmPixelHit* hit) const { return GetTrackingStationIndex(hit->GetAddress()); }
+
+ /// Gets a tracking station by the address
+ /// \param address Unique element address
+ /// \return Local index of the tracking station
+ int GetTrackingStationIndex(int address) const
+ {
+ return CbmTrdAddress::GetLayerId(address);
+ }
/// Gets max size of a station along the X-axis
/// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
diff --git a/core/qa/CbmQaHist.h b/core/qa/CbmQaHist.h
index 860e9d1e8b..edb6f195ce 100644
--- a/core/qa/CbmQaHist.h
+++ b/core/qa/CbmQaHist.h
@@ -84,7 +84,7 @@ public:
f->SetLineColor(kRed);
f->SetLineWidth(3);
TPaveStats* st = (TPaveStats*) this->FindObject("stats");
- if (!st) { LOG(fatal) << "CbmQaHist: can not access histogram statistics"; }
+ if (!st) { LOG(fatal) << "CbmQaHist: can not access statistics of histogram with name \"" << this->GetName() << '\"'; }
else {
st->SetX1NDC(0.6);
st->SetX2NDC(0.940);
diff --git a/mvd/CbmMvdTrackingInterface.h b/mvd/CbmMvdTrackingInterface.h
index b2df5b8418..01b845e3ed 100644
--- a/mvd/CbmMvdTrackingInterface.h
+++ b/mvd/CbmMvdTrackingInterface.h
@@ -12,6 +12,7 @@
#ifndef CbmMvdTrackingInterface_h
#define CbmMvdTrackingInterface_h 1
+#include "CbmMvdDetectorId.h"
#include "CbmMvdDetector.h"
#include "CbmMvdHit.h"
#include "CbmMvdStationPar.h"
@@ -28,7 +29,7 @@
/// Class CbmMvdTrackingInterface is a CbmL1 subtask, which provides necessary methods for L1 tracker
/// to access the geometry and dataflow settings.
///
-class CbmMvdTrackingInterface : public FairTask, public CbmTrackingDetectorInterfaceBase {
+class CbmMvdTrackingInterface : public FairTask, public CbmTrackingDetectorInterfaceBase, public CbmMvdDetectorId {
public:
/// Default constructor
CbmMvdTrackingInterface();
@@ -115,6 +116,14 @@ public:
}();
return hitMvd->GetStationNr();
}
+
+ /// Gets a tracking station by the address of element (detectorID in terms of MVD)
+ /// \param detectorId Unique element address (detectorID in terms of MVD)
+ /// \return Local index of the tracking station
+ int GetTrackingStationIndex(int detectorId) const
+ {
+ return StationNr(detectorId);
+ }
/// Gets max size of a tracking station along the X-axis
/// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
diff --git a/reco/L1/CMakeLists.txt b/reco/L1/CMakeLists.txt
index f21620c1f0..4f43b55df7 100644
--- a/reco/L1/CMakeLists.txt
+++ b/reco/L1/CMakeLists.txt
@@ -157,7 +157,6 @@ ParticleFinder/CbmL1PFFitter.cxx
ParticleFinder/CbmL1PFMCParticle.cxx
qa/CbmTrackerInputQaTrd.cxx
-qa/CbmTrackingInputQaBase.cxx
qa/CbmTrackingInputQaSts.cxx
)
@@ -191,7 +190,6 @@ OffLineInterface/CbmL1GlobalFindTracksEvents.h
L1Algo/L1Def.h
L1Algo/L1Vector.h
qa/CbmTrackerInputQaTrd.h
-qa/CbmTrackingInputQaBase.h
qa/CbmTrackingInputQaSts.h
)
diff --git a/reco/L1/CbmL1ReadEvent.cxx b/reco/L1/CbmL1ReadEvent.cxx
index e74aedc781..45ea39e4cf 100644
--- a/reco/L1/CbmL1ReadEvent.cxx
+++ b/reco/L1/CbmL1ReadEvent.cxx
@@ -1482,8 +1482,8 @@ void CbmL1::HitMatch()
static_cast<const CbmMatch*>(listStsClusterMatch->At(sh->GetBackClusterId()));
CbmMatch stsHitMatch;
- Float_t Sum_of_front = 0;
- Float_t Sum_of_back = 0;
+ Float_t Sum_of_front = 0; // total weight of all front links
+ Float_t Sum_of_back = 0; // total weight of all back links
for (Int_t iFrontLink = 0; iFrontLink < frontClusterMatch->GetNofLinks(); iFrontLink++) {
diff --git a/reco/L1/L1LinkDef.h b/reco/L1/L1LinkDef.h
index a4407477ec..7bb7c97cbf 100644
--- a/reco/L1/L1LinkDef.h
+++ b/reco/L1/L1LinkDef.h
@@ -29,7 +29,6 @@
//#pragma link C++ class CbmL1SttHit+;
//#pragma link C++ class CbmL1SttTrackFinder+;
//#pragma link C++ class CbmL1SttTrack+;
-#pragma link C++ class CbmTrackingInputQaBase + ;
#pragma link C++ class CbmTrackingInputQaSts + ;
#pragma link C++ class CbmTrackerInputQaTrd + ;
diff --git a/reco/L1/qa/CbmTrackingInputQaBase.cxx b/reco/L1/qa/CbmTrackingInputQaBase.cxx
deleted file mode 100644
index d0ffaf2475..0000000000
--- a/reco/L1/qa/CbmTrackingInputQaBase.cxx
+++ /dev/null
@@ -1,319 +0,0 @@
-/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
- SPDX-License-Identifier: GPL-3.0-only
- Authors: Sergey Gorbunov, Sergei Zharko [committer] */
-
-/***************************************************************************************************
- * @file CbmTrackingInputQaBase.cxx
- * @brief Base class for the tracking input QA (definition)
- * @since 24.06.2022
- * @author S.Zharko <s.zharko@gsi.de>
- ***************************************************************************************************/
-
-#include "CbmTrackingInputQaBase.h"
-
-#include "CbmCluster.h"
-#include "CbmDigiManager.h"
-#include "CbmMCEventList.h"
-#include "CbmMatch.h"
-#include "CbmPixelHit.h"
-#include "CbmStsAddress.h" // TODO: TMP for tests, must be removed!!!! (S.Zharko)
-#include "CbmTimeSlice.h"
-
-#include "FairLogger.h"
-
-#include "TClonesArray.h"
-
-#include <iostream>
-
-// --------------------------------------------------------------------------------------------------------------------
-//
-CbmTrackingInputQaBase::CbmTrackingInputQaBase(const char* name, int verbose) : FairTask(name, verbose) {}
-
-// --------------------------------------------------------------------------------------------------------------------
-//
-CbmTrackingInputQaBase::~CbmTrackingInputQaBase() {}
-
-ClassImp(CbmTrackingInputQaBase);
-
-//
-// --------------------------------------------------------------------------------------------------------------------
-//
-void CbmTrackingInputQaBase::Clear(Option_t*)
-{
- std::cout << "\033[1;32mCALL:\033[0mCbmTrackingInputQaBase::Clear(Option_t*)\n";
-}
-//
-// --------------------------------------------------------------------------------------------------------------------
-//
-void CbmTrackingInputQaBase::Exec(Option_t*)
-{
- std::cout << "\033[1;32mCALL:\033[0mCbmTrackingInputQaBase::Exec(Option_t*)\n";
-
- // Run resolution Qa
- ResolutionQa();
-
- // Update number of events
- fNevents.SetVal(fNevents.GetVal() + 1);
-}
-//
-// --------------------------------------------------------------------------------------------------------------------
-//
-void CbmTrackingInputQaBase::Finish() { std::cout << "\033[1;32mCALL:\033[0mCbmTrackingInputQaBase::Finish()\n"; }
-//
-// --------------------------------------------------------------------------------------------------------------------
-//
-InitStatus CbmTrackingInputQaBase::Init()
-{
- std::cout << "\033[1;32mCALL:\033[0mCbmTrackingInputQaBase::Init()\n";
-
- // If detector is not present, fpDetectorInterface is nullptr
- if (!fpDetectorInterface) { LOG(fatal) << "CbmTrackingInputQa: Tracking detector interface is undefined"; }
-
- // Input data initialization
- InitInputData();
-
- // Base histograms initialization
- InitBaseHistograms();
-
- // Register histograms in output file
- // TODO: Add subfolders for each subsystem like
- fpOutFolder =
- std::make_shared<TFolder>(TString("TrackingInputQa_") + fDetName, TString("TrackingInputQa_") + fDetName);
- fpOutFolder->SetOwner(false);
- fpOutFolderHists = fpOutFolder->AddFolder("rawHist", "Raw histograms");
- gStyle->SetOptStat(0);
-
- fNevents.SetVal(0);
- fpOutFolderHists->Add(&fNevents);
-
- // Register histograms in the output folder
- for (auto* histo : fHistList) {
- fpOutFolderHists->Add(histo);
- }
-
- // Geometry QA
- GeometryQa();
-
- return kSUCCESS;
-}
-
-// --------------------------------------------------------------------------------------------------------------------
-//
-InitStatus CbmTrackingInputQaBase::InitInputData()
-{
- std::cout << "\033[1;32mCALL:\033[0mCbmTrackingInputQaBase::InitInputData()\n";
-
- auto* manager = FairRootManager::Instance();
- fpDigiManager = CbmDigiManager::Instance();
- fpDigiManager->Init(); // NOTE: is initialized only once
-
- // Reconstructed data initialization
- fpTimeSlice = dynamic_cast<CbmTimeSlice*>(manager->GetObject("TimeSlice."));
- if (!fpTimeSlice) {
- LOG(error) << "Time slice was not found";
- return kERROR;
- }
-
- fpHits = GetHitsInput();
- if (!fpHits) {
- LOG(error) << "Hits input array was not found for " << fDetTitle;
- return kERROR;
- }
-
- fpClusters = GetClustersInput();
- if (!fpClusters) {
- LOG(error) << "Cluster input array was not found for " << fDetTitle;
- return kERROR;
- }
-
- // MC data initialization
- fpMcManager = dynamic_cast<CbmMCDataManager*>(manager->GetObject("MCDataManager"));
- fIsMcPresent = bool(fpMcManager);
-
- if (fIsMcPresent) {
- fpMcEventList = dynamic_cast<CbmMCEventList*>(manager->GetObject("MCEventList."));
- fpMcTracks = fpMcManager->InitBranch("MCTrack");
- fpMcPoints = GetMcPointsInput();
- fpHitMatches = GetHitMatchesInput();
-
- if (!fpMcEventList) {
- LOG(error) << "MCEventList data were not found";
- return kERROR;
- }
-
- if (!fpMcTracks) {
- LOG(error) << "MC tracks data were not found";
- return kERROR;
- }
-
- if (!fpMcPoints) {
- LOG(error) << "MC points data were not found for " << fDetTitle;
- return kERROR;
- }
-
- if (!fpHitMatches) {
- LOG(error) << "Hits match data were not found for " << fDetTitle;
- return kERROR;
- }
-
- if (!fpDigiManager->IsMatchPresent(fDetId)) {
- LOG(error) << "No digi matches were found for " << fDetTitle;
- return kERROR;
- }
- }
-
- return kSUCCESS;
-}
-
-// --------------------------------------------------------------------------------------------------------------------
-//
-InitStatus CbmTrackingInputQaBase::InitBaseHistograms()
-{
- std::cout << "\033[1;32mCALL:\033[0mCbmTrackingInputQaBase::InitBaseHistograms()\n";
-
- fHistResidualX.clear();
- fHistResidualY.clear();
- fHistResidualT.clear();
-
- fHistPullX.clear();
- fHistPullY.clear();
- fHistPullT.clear();
-
- fHistPointsPerHit.clear();
- fHistHitsPerPoint.clear();
-
- int nTrackingStations = fpDetectorInterface->GetNtrackingStations();
-
- fHistResidualX.reserve(nTrackingStations);
- fHistResidualY.reserve(nTrackingStations);
- fHistResidualT.reserve(nTrackingStations);
-
- fHistPullX.reserve(nTrackingStations);
- fHistPullY.reserve(nTrackingStations);
- fHistPullT.reserve(nTrackingStations);
-
- fHistPointsPerHit.reserve(nTrackingStations);
- fHistHitsPerPoint.reserve(nTrackingStations);
-
- // Create histograms
- TString histName = "";
- TString histTitle = "";
- int nBins = -1;
- int rMin = 0.;
- int rMax = 0.;
- for (int iSt = 0; iSt < nTrackingStations; ++iSt) {
- // Residuals
- histName = fDetName + Form("_st%d_h1ResidualX", iSt);
- histTitle = fDetTitle + Form(" Station %d: Residual X; (X_{reco} - X_{MC}) [cm]", iSt);
- std::tie(nBins, rMin, rMax) = fRangeResidualX;
- fHistResidualX.emplace_back(histName, histTitle, nBins, rMin, rMax);
-
- histName = fDetName + Form("_st%d_h1ResidualY", iSt);
- histTitle = fDetTitle + Form(" Station %d: Residual Y; (Y_{reco} - Y_{MC}) [cm]", iSt);
- std::tie(nBins, rMin, rMax) = fRangeResidualY;
- fHistResidualY.emplace_back(histName, histTitle, nBins, rMin, rMax);
-
- histName = fDetName + Form("_st%d_h1ResidualT", iSt);
- histTitle = fDetTitle + Form(" Station %d: Residual T; (T_{reco} - T_{MC}) [cm]", iSt);
- std::tie(nBins, rMin, rMax) = fRangeResidualT;
- fHistResidualT.emplace_back(histName, histTitle, nBins, rMin, rMax);
-
- // Pulls
- histName = fDetName + Form("_st%d_h1PullX", iSt);
- histTitle = fDetTitle + Form(" Station %d: Pull X; (X_{reco} - X_{MC}) / #sigmaX_{reco}", iSt);
- std::tie(nBins, rMin, rMax) = fRangePullX;
- fHistPullX.emplace_back(histName, histTitle, nBins, rMin, rMax);
-
- histName = fDetName + Form("_st%d_h1PullY", iSt);
- histTitle = fDetTitle + Form(" Station %d: Pull Y; (Y_{reco} - Y_{MC}) / #sigmaY_{reco}", iSt);
- std::tie(nBins, rMin, rMax) = fRangePullY;
- fHistPullY.emplace_back(histName, histTitle, nBins, rMin, rMax);
-
- histName = fDetName + Form("_st%d_h1PullT", iSt);
- histTitle = fDetTitle + Form(" Station %d: Pull T; (T_{reco} - T_{MC}) / #sigmaT_{reco}", iSt);
- std::tie(nBins, rMin, rMax) = fRangePullT;
- fHistPullT.emplace_back(histName, histTitle, nBins, rMin, rMax);
-
- // MC points per one hit
- histName = fDetName + Form("_st%d_h1PointsPerHit", iSt);
- histTitle = fDetTitle + Form(" Station %d: MC Points per Hit; N_{MC Points} per hit", iSt);
- std::tie(nBins, rMin, rMax) = fRangePointsPerHit;
- fHistPointsPerHit.emplace_back(histName, histTitle, nBins, rMin, rMax);
-
- // Hits per one mc point
- histName = fDetName + Form("_st%d_h1HitsPerPoint", iSt);
- histTitle = fDetTitle + Form(" Station %d: Hit per MC points; N_{hits} per MC point", iSt);
- std::tie(nBins, rMin, rMax) = fRangeHitsPerPoint;
- fHistHitsPerPoint.emplace_back(histName, histTitle, nBins, rMin, rMax);
- }
-
- // Register histograms
- for (int iSt = 0; iSt < nTrackingStations; ++iSt) {
- RegisterHist(&fHistResidualX[iSt]);
- RegisterHist(&fHistResidualY[iSt]);
- RegisterHist(&fHistResidualT[iSt]);
-
- RegisterHist(&fHistPullX[iSt]);
- RegisterHist(&fHistPullY[iSt]);
- RegisterHist(&fHistPullT[iSt]);
-
- RegisterHist(&fHistPointsPerHit[iSt]);
- RegisterHist(&fHistHitsPerPoint[iSt]);
- }
-
- return kSUCCESS;
-}
-
-// --------------------------------------------------------------------------------------------------------------------
-//
-void CbmTrackingInputQaBase::Reset() { std::cout << "\033[1;32mCALL:\033[0mCbmTrackingInputQaBase::Reset()\n"; }
-
-// --------------------------------------------------------------------------------------------------------------------
-//
-void CbmTrackingInputQaBase::SetParContainers()
-{
- std::cout << "\033[1;32mCALL:\033[0mCbmTrackingInputQaBase::SetParContainers()\n";
-}
-
-// --------------------------------------------------------------------------------------------------------------------
-//
-InitStatus CbmTrackingInputQaBase::GeometryQa()
-{
- std::cout << "\033[1;32mCALL:\033[0mCbmTrackingInputQaBase::GeometryQa()\n";
- return kSUCCESS;
-}
-
-// --------------------------------------------------------------------------------------------------------------------
-//
-void CbmTrackingInputQaBase::ResolutionQa()
-{
- std::cout << "\033[1;32mCALL:\033[0mCbmTrackingInputQaBase::ResolutionQa()\n";
-
- // Resolution QA is impossible without MC information
- if (!fIsMcPresent) { return; }
-
- int nHits = fpHits->GetEntriesFast(); // Hits in the event/time slice
- int nClusters = fpClusters->GetEntriesFast(); // Clusters in the event/time slice
- int nDigis = fpDigiManager->GetNofDigis(fDetId); // Number of digies
- int nStations = fpDetectorInterface->GetNtrackingStations(); // Number of tracking stations
-
- // *********************************************
- // ** Loop over hits in this event/time slice **
- // *********************************************
-
- for (int iHit = 0; iHit < nHits; ++iHit) {
- /// TODO: put it into a function from iHit
-
- auto* hit = dynamic_cast<CbmPixelHit*>(fpHits->At(iHit));
- if (!hit) { LOG(error) << "There is an empty hit at position " << iHit << " (" << fDetTitle << ")"; }
-
- // Tracking station index
- int iStation = fpDetectorInterface->GetTrackingStationIndex(hit);
- if (iStation < 0 || iStation >= nStations) {
- LOG(fatal) << "Wrong tracking station index was calculated by the tracking detector interface for " << fDetTitle
- << ": station index = " << iStation << " (expected in the range [0, " << nStations - 1 << "])";
- return;
- }
-
- } // Loop over hits in this event/time slice: END
-}
diff --git a/reco/L1/qa/CbmTrackingInputQaBase.h b/reco/L1/qa/CbmTrackingInputQaBase.h
deleted file mode 100644
index 76feaf7ba9..0000000000
--- a/reco/L1/qa/CbmTrackingInputQaBase.h
+++ /dev/null
@@ -1,195 +0,0 @@
-/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
- SPDX-License-Identifier: GPL-3.0-only
- Authors: Sergey Gorbunov, Sergei Zharko [committer] */
-
-/***************************************************************************************************
- * @file CbmTrackingInputQaBase.h
- * @brief Base class for the tracking input QA (declaration)
- * @since 24.06.2022
- * @author S.Zharko <s.zharko@gsi.de>
- ***************************************************************************************************/
-
-#ifndef CbmTrackingInputQaBase_h
-#define CbmTrackingInputQaBase_h 1
-
-#include "CbmMCDataManager.h"
-#include "CbmQaCanvas.h"
-#include "CbmQaHist.h"
-#include "CbmSetup.h"
-#include "CbmTrackingDetectorInterfaceBase.h" // Communicate via tracking detector interface
-
-#include <FairTask.h>
-
-#include "TFolder.h"
-#include "TH1D.h"
-#include "TH1I.h"
-#include "TH2D.h"
-#include "TParameter.h"
-#include "TProfile.h"
-#include "TString.h"
-#include <tuple>
-
-#include <memory>
-#include <vector>
-
-class CbmDigiManager;
-class CbmMCEventList;
-class CbmMCDataArray;
-class CbmTimeSlice;
-class TClonesArray;
-class CbmCluster;
-class CbmMatch;
-
-/// Class CbmTrackingInputQaBase is an abstract FairTask inherited class, providing interface for a
-/// tracking detector QA including input geometry and data checks. One should inherit this class
-/// for a particular detector subsystem.
-///
-class CbmTrackingInputQaBase : public FairTask {
-public:
- /// Constructor from the task name and verbosity level
- CbmTrackingInputQaBase(const char* name = "CbmTrackingInputQaBase", int verbose = 1);
-
- /// Destructor
- virtual ~CbmTrackingInputQaBase();
-
- /// Suppress copy and move
- CbmTrackingInputQaBase(const CbmTrackingInputQaBase&) = delete;
- CbmTrackingInputQaBase(CbmTrackingInputQaBase&&) = delete;
- CbmTrackingInputQaBase& operator=(const CbmTrackingInputQaBase&) = delete;
- CbmTrackingInputQaBase& operator=(CbmTrackingInputQaBase&&) = delete;
-
- /// TTask: Clears all data structures created by a previous execution of the task
- /// \param option Specify action for the TTask::Exec (not defined yet)
- void Clear(Option_t* option = "");
-
- /// Ttask: Executes the task (processes a timeslice)
- /// \param option Specify action for the TTask::Exec (not defined yet)
- void Exec(Option_t* /*option*/);
-
- /// FairTask: Action at hte end of the run
- void Finish();
-
- //TFolder& GetQa();
-
- /// FairTask: Task initialization an the beginning of the run
- InitStatus Init();
-
- /// FairTask: Task reinitialization
- InitStatus ReInit() { return Init(); }
-
- /// FairTask: Set parameter containers (if any)
- void SetParContainers();
-
- /// Resets data fields
- void Reset();
-
- /// Gets pointer to hits array of a particular detector type
- virtual TClonesArray* GetHitsInput() = 0;
-
- /// Gets pointer to clusters array of a particular detector type
- virtual TClonesArray* GetClustersInput() = 0;
-
- /// Gets pointer to hit matches array for a particular detector type
- virtual TClonesArray* GetHitMatchesInput() = 0;
-
- /// Gets pointer to MC points for a particular detector type
- virtual CbmMCDataArray* GetMcPointsInput() = 0;
-
- ClassDef(CbmTrackingInputQaBase, 0); // Base class for tracking input QA
-
-protected:
- /// Registers histogram
- void RegisterHist(TH1* histo)
- {
- fHistList.push_back(histo);
- histo->SetDirectory(0);
- }
-
- /// Checks geometry
- InitStatus GeometryQa(); // TODO: Final or not final? (S.Zharko)
-
- /// Checks resolution
- /// A general scheme of resolution checks is provided in this method. One should specify particular operations in the
- /// derived classes for each particular detector....
- void ResolutionQa(); // TODO: Final or not final? (S.Zharko)
-
-
- // ***********
- // ** Setup **
- // ***********
-
-
- CbmTrackingDetectorInterfaceBase* fpDetectorInterface {nullptr}; ///< Pointer to current tracking detector I/F
- CbmTimeSlice* fpTimeSlice {nullptr}; ///< Pointer to current time slice
- CbmDigiManager* fpDigiManager {nullptr}; ///< Pointer to digi manager
-
- CbmMCEventList* fpMcEventList {nullptr};
- CbmMCDataManager* fpMcManager {nullptr};
-
- // **********
- // ** Data **
- // **********
-
- TClonesArray* fpClusters {nullptr};
- TClonesArray* fpHits {nullptr};
- TClonesArray* fpHitMatches {nullptr};
-
- CbmMCDataArray* fpMcTracks {nullptr};
- CbmMCDataArray* fpMcPoints {nullptr};
-
- // ********************************
- // ** List of base QA histograms **
- // ********************************
-
- // NOTE: One can add extra detector specific histograms in the Init method of the derived class. In that case, the
- // histograms should be registered via RegisterHisto method
-
- std::vector<CbmQaHist<TH1D>> fHistResidualX; ///< Residual distributions for X vs index of a station
- std::vector<CbmQaHist<TH1D>> fHistResidualY; ///< Residual distributions for X vs index of a station
- std::vector<CbmQaHist<TH1D>> fHistResidualT; ///< Residual distributions for T vs index of a station
-
- std::vector<CbmQaHist<TH1D>> fHistPullX; ///< Pulls distribution for X vs. index of a station
- std::vector<CbmQaHist<TH1D>> fHistPullY; ///< Pulls distribution for Y vs. index of a station
- std::vector<CbmQaHist<TH1D>> fHistPullT; ///< Pulls distribution for T vs. index of a station
-
- std::vector<CbmQaHist<TH1I>> fHistPointsPerHit; ///< Distribution of MC points per one hit vs. index of a station
- std::vector<CbmQaHist<TH1I>> fHistHitsPerPoint; ///< Distribution of hits per one MC point vs. index of a station
-
- // Histogram properties, can be modified in the constructor of a derived class
- std::tuple<int, double, double> fRangeResidualX {100, -5., 5.};
- std::tuple<int, double, double> fRangeResidualY {100, -10., 10.};
- std::tuple<int, double, double> fRangeResidualT {100, -50., 50.};
-
- std::tuple<int, double, double> fRangePullX {100, -5., 5.};
- std::tuple<int, double, double> fRangePullY {100, -5., 5.};
- std::tuple<int, double, double> fRangePullT {100, -5., 5.};
-
- std::tuple<int, double, double> fRangePointsPerHit {10, -0.5, 9.5};
- std::tuple<int, double, double> fRangeHitsPerPoint {10, -0.5, 9.5};
-
- TParameter<long> fNevents {"nEvents", 0}; ///< Number of processed events
-
- // ************
- // ** Output **
- // ************
-
- std::shared_ptr<TFolder> fpOutFolder {nullptr}; ///< Output folder for this QA task
- TFolder* fpOutFolderHists {nullptr}; ///< Subfolder for raw histograms
-
- ECbmModuleId fDetId; ///< ID of the detector subsystem
- TString fDetName {""}; ///< Name of the detector subsystem, used for i/o (e.g., mvd, sts, much, trd1d, trd2d, tof)
- TString fDetTitle {""}; ///< Title of the detector subsystem, used for logs (e.g., MVD, STS, MuCh, TRD-1D, ...)
-
-private:
- /// Initializes base histograms
- InitStatus InitBaseHistograms();
-
- /// Initializes input data arrays
- InitStatus InitInputData();
-
- std::vector<TH1*> fHistList; ///< List of the pointers to all histograms contained in the class
-
- bool fIsMcPresent {false};
-};
-
-#endif // CbmTrackingInputQaBase_h
diff --git a/reco/L1/qa/CbmTrackingInputQaSts.cxx b/reco/L1/qa/CbmTrackingInputQaSts.cxx
index b02229e00f..29a546ef9d 100644
--- a/reco/L1/qa/CbmTrackingInputQaSts.cxx
+++ b/reco/L1/qa/CbmTrackingInputQaSts.cxx
@@ -4,66 +4,799 @@
/***************************************************************************************************
* @file CbmTrackingInputQaSts.cxx
- * @brief Class for the tracking input QA for STS (definition)
- * @since 27.06.2022
+ * @brief Base class for the tracking input QA (definition)
+ * @since 24.06.2022
* @author S.Zharko <s.zharko@gsi.de>
***************************************************************************************************/
#include "CbmTrackingInputQaSts.h"
-#include "CbmDigiManager.h"
#include "CbmMCDataArray.h"
+#include "CbmStsCluster.h"
+#include "CbmDigiManager.h"
#include "CbmMCEventList.h"
#include "CbmMatch.h"
-#include "CbmSetup.h"
-#include "CbmStsCluster.h"
-#include "CbmStsDigi.h"
#include "CbmStsHit.h"
-#include "CbmStsTrackingInterface.h"
+#include "CbmStsAddress.h" // TODO: TMP for tests, must be removed!!!! (S.Zharko)
+#include "CbmTimeSlice.h"
+#include "CbmMCTrack.h"
+#include "CbmStsPoint.h"
+#include "CbmStsTrackingInterface.h" // Communicate via tracking detector interface
+
+#include "FairLogger.h"
+
+#include "TClonesArray.h"
+#include "TParticlePDG.h"
+#include "TDatabasePDG.h"
+
+#include <iostream>
+#include <stdexcept>
ClassImp(CbmTrackingInputQaSts);
-// --------------------------------------------------------------------------------------------------------------------
+// ---------------------------------------------------------------------------------------------------------------------
//
-CbmTrackingInputQaSts::CbmTrackingInputQaSts(int verbose) : CbmTrackingInputQaBase("CbmTrackingInputQaSts", verbose)
+CbmTrackingInputQaSts::CbmTrackingInputQaSts(int verbosity) : FairTask("CbmTrackingInputQaSts", verbosity)
{
- std::cout << "\033[1;32mCALL:\033[0mCbmTrackingInputQaSts::CbmTrackingInputQaSts()\n";
+}
- // Setup detector-specific parameters
- fDetId = ECbmModuleId::kSts;
- fDetName = "sts";
- fDetTitle = "STS";
- fpDetectorInterface = CbmStsTrackingInterface::Instance();
+// ---------------------------------------------------------------------------------------------------------------------
+//
+CbmTrackingInputQaSts::~CbmTrackingInputQaSts()
+{
+ DeInit();
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+bool CbmTrackingInputQaSts::CheckDistributions()
+{
+ bool res = true;
+
+ const int nStations = fpDetectorInterface->GetNtrackingStations();
+
+ // ** Check pulls **
+ //
+
+
+ // Function to compare pulls distributions RMS with 1
+ // TODO: Choose proper selection criteria (S.Zharko)
+ auto CheckPullsDistribution = [&](const TH1& pHist) {
+ // Checks sigma of distribution fit with unity
+ auto* pFitFunc = pHist.GetFunction("gaus");
+ if (!pFitFunc) {
+ throw std::runtime_error("STS tracking input QA: attempt to check sigma of unfitted pulls histogram");
+ }
+ auto vSigma = pFitFunc->GetParameter(2);
+ auto eSigma = pFitFunc->GetParError(2);
+
+ // Select 3 sigma interval
+ LOG(info) << "Checking histogram \"" << pHist.GetName() << "\": fit result = " << vSigma << " +/- " << eSigma;
+ if (std::fabs(vSigma - 1.) < 3. * eSigma) {
+ return true;
+ }
+ else {
+ return false;
+ }
+
+ //auto vRms = pHist.GetRMS(); // RMS value
+ //auto eRms = pHist.GetRMSError(); // RMS error
+ //std::cout << std::fabs(vRms - 1.) << '\n';
+ //if (eRms > vRms * fMaxPullsRmsDiff ) {
+ // LOG(error) << "STS tracking input QA: the analysis of \"" << pHist.GetName() << "\" RMS is impossible due to "
+ // << "lack of statistics\n";
+ // // TODO: probably, here we can use fit? (S.Zharko)
+ // return false;
+ //}
+ //if (std::fabs(vRms - 1.) < fMaxPullsRmsDiff) {
+ // return true;
+ //}
+ //else {
+ // LOG(info) << "STS tracking input QA: RMS check for pulls distribution failed (histogram: \"" << pHist.GetName() << "\"): "
+ // << "RMS = " << vRms << " +/- " << eRms;
+ // return false;
+ //}
+ return false;
+ };
+
+
+ for (int iSt = 0; iSt < nStations; ++iSt) {
+ res = CheckPullsDistribution(fHistPullX[iSt]) && res;
+ res = CheckPullsDistribution(fHistPullY[iSt]) && res;
+ res = CheckPullsDistribution(fHistPullT[iSt]) && res;
+ }
+
+ return res;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void CbmTrackingInputQaSts::DeInit()
+{
+ fIsMcPresent = false;
+
+ fpDetectorInterface = nullptr;
+ fpTimeSlice = nullptr;
+ fpDigiManager = nullptr;
+
+ fpMcManager = nullptr;
+ fpMcTracks = nullptr;
+ fpMcPoints = nullptr;
+
+ fpMcEventList = nullptr;
+ fpClusters = nullptr;
+ fpHits = nullptr;
+ fpClusterMatches = nullptr;
+ fpHitMatches = nullptr;
+
+ fOutFolder.Clear();
+ fpOutFolderHists = nullptr;
+
+ fNevents.SetVal(0);
+
+ fHistResidualX.clear();
+ fHistResidualY.clear();
+ fHistResidualT.clear();
+
+ fHistPullX.clear();
+ fHistPullY.clear();
+ fHistPullT.clear();
+
+ fHistPointsPerHit.clear();
+ fHistHitsPerPoint.clear();
+
+ fHistEfficiencyXY.clear();
+ fHistEfficiencyR.clear();
}
// --------------------------------------------------------------------------------------------------------------------
//
-CbmTrackingInputQaSts::~CbmTrackingInputQaSts() {}
+void CbmTrackingInputQaSts::Exec(Option_t*)
+{
+ // Run resolution Qa
+ ResolutionQa();
+
+ // Update number of events
+ fNevents.SetVal(fNevents.GetVal() + 1);
+}
// --------------------------------------------------------------------------------------------------------------------
//
+void CbmTrackingInputQaSts::Finish()
+{
+ // Add output to a sink
+ auto* pSink = FairRootManager::Instance()->GetSink();
+ if (pSink) {
+ pSink->WriteObject(&GetQa(), nullptr);
+ }
+
+ // Check accumulated distributions
+ bool areResolutionsOk = CheckDistributions();
+
+ // TODO: Collect all the flags in one place and make a decission here (S.Zharko)
+
+ if (areResolutionsOk) {
+ LOG(info) << this->GetName() << ": \033[1;32mtask succeeded\033[0m";
+ }
+ else {
+ LOG(info) << this->GetName() << ": \033[1;31mtask failed\033[0m";
+ }
+
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+InitStatus CbmTrackingInputQaSts::GeometryQa()
+{
+ // Test stations ordering: TODO - is it needed here? This test is done in L1
+ return kSUCCESS;
+}
+
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+TFolder& CbmTrackingInputQaSts::GetQa()
+{
+ //gStyle->SetPaperSize(20, 20);
+
+ // Loop over tracking stations
+ const int nStations = fpDetectorInterface->GetNtrackingStations();
+ for (int iSt = 0; iSt < nStations; ++iSt) {
+ // Fit histograms
+ fHistResidualX[iSt].Fit("gaus", "Q");
+ fHistResidualY[iSt].Fit("gaus", "Q");
+ fHistResidualT[iSt].Fit("gaus", "Q");
+
+ fHistPullX[iSt].Fit("gaus", "Q");
+ fHistPullY[iSt].Fit("gaus", "Q");
+ fHistPullT[iSt].Fit("gaus", "Q");
+
+ // Draw histograms
+ fCanvResidualX.cd(iSt + 1);
+ fHistResidualX[iSt].DrawCopy("", "");
+ fCanvResidualY.cd(iSt + 1);
+ fHistResidualY[iSt].DrawCopy("", "");
+ fCanvResidualT.cd(iSt + 1);
+ fHistResidualT[iSt].DrawCopy("", "");
+
+ fCanvPullX.cd(iSt + 1);
+ fHistPullX[iSt].DrawCopy("", "");
+ fCanvPullY.cd(iSt + 1);
+ fHistPullY[iSt].DrawCopy("", "");
+ fCanvPullT.cd(iSt + 1);
+ fHistPullT[iSt].DrawCopy("", "");
+
+ fCanvPointsPerHit.cd(iSt + 1);
+ fHistPointsPerHit[iSt].DrawCopy("", "");
+
+ fCanvHitsPerPoint.cd(iSt + 1);
+ fHistHitsPerPoint[iSt].DrawCopy("", "");
+
+ fCanvEfficiencyR.cd(iSt + 1);
+ fHistEfficiencyR[iSt].DrawCopy("", "");
+
+ fCanvEfficiencyXY.cd(iSt + 1);
+ fHistEfficiencyXY[iSt].DrawCopy("colz", "");
+ }// Loop over tracking stations: end
+ return fOutFolder;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
InitStatus CbmTrackingInputQaSts::Init()
{
- std::cout << "\033[1;32mCALL:\033[0mCbmTrackingInputQaSts::Init()\n";
- // Check, if STS is in the setup
- if (!CbmSetup::Instance()->IsActive(ECbmModuleId::kSts)) {
- LOG(warn) << "STS is not active in the current setup";
- return kSUCCESS;
+ // Reset data fields
+ DeInit();
+
+ // If detector is not present, fpDetectorInterface is nullptr
+ fpDetectorInterface = CbmStsTrackingInterface::Instance();
+ if (!fpDetectorInterface) { LOG(fatal) << "CbmTrackingInputQa: Tracking detector interface is undefined"; }
+
+ // Input data initialization
+ ReadAndCreateDataBranches();
+
+ // Histograms initialization
+ InitHistograms();
+
+ // Canvases initialization
+ InitCanvases();
+
+ // Register histograms in output file
+ // TODO: Add subfolders for each subsystem like
+ fOutFolder.SetOwner(false);
+ fpOutFolderHists = fOutFolder.AddFolder("rawHist", "Raw histograms");
+ gStyle->SetOptStat(0);
+
+ fNevents.SetVal(0); // redundant
+ fpOutFolderHists->Add(&fNevents);
+
+ // Register histograms in the output folder
+ for (auto* histo : fHistList) {
+ fpOutFolderHists->Add(histo);
}
+
+ // Do QA of geometry in the end of initialization
+ return GeometryQa();
+}
+
+
+// --------------------------------------------------------------------------------------------------------------------
+//
+InitStatus CbmTrackingInputQaSts::InitCanvases()
+{
+ // Reset canvases
+ fCanvResidualX.Clear();
+ fCanvResidualY.Clear();
+ fCanvResidualT.Clear();
+
+ fCanvPullX.Clear();
+ fCanvPullY.Clear();
+ fCanvPullT.Clear();
+
+ fCanvPointsPerHit.Clear();
+ fCanvHitsPerPoint.Clear();
+
+ fCanvEfficiencyR.Clear();
+ fCanvEfficiencyXY.Clear();
- // Initialize base histograms
- CbmTrackingInputQaBase::Init();
+ // Devide canvases into sections to store plots vs. station ID
+ fCanvResidualX.Divide2D(fpDetectorInterface->GetNtrackingStations());
+ fCanvResidualY.Divide2D(fpDetectorInterface->GetNtrackingStations());
+ fCanvResidualT.Divide2D(fpDetectorInterface->GetNtrackingStations());
- // Additional STS-specific histograms can be initialized here
- // ...
+ fCanvPullX.Divide2D(fpDetectorInterface->GetNtrackingStations());
+ fCanvPullY.Divide2D(fpDetectorInterface->GetNtrackingStations());
+ fCanvPullT.Divide2D(fpDetectorInterface->GetNtrackingStations());
+ fCanvPointsPerHit.Divide2D(fpDetectorInterface->GetNtrackingStations());
+ fCanvHitsPerPoint.Divide2D(fpDetectorInterface->GetNtrackingStations());
+
+ fCanvEfficiencyR.Divide2D(fpDetectorInterface->GetNtrackingStations());
+ fCanvEfficiencyXY.Divide2D(fpDetectorInterface->GetNtrackingStations());
+
+ // Add the canvases to the output folder
+ fOutFolder.Add(&fCanvResidualX);
+ fOutFolder.Add(&fCanvResidualY);
+ fOutFolder.Add(&fCanvResidualT);
+
+ fOutFolder.Add(&fCanvPullX);
+ fOutFolder.Add(&fCanvPullY);
+ fOutFolder.Add(&fCanvPullT);
+
+ fOutFolder.Add(&fCanvPointsPerHit);
+ fOutFolder.Add(&fCanvHitsPerPoint);
+
+ fOutFolder.Add(&fCanvEfficiencyR);
+ fOutFolder.Add(&fCanvEfficiencyXY);
+
return kSUCCESS;
}
// --------------------------------------------------------------------------------------------------------------------
//
+InitStatus CbmTrackingInputQaSts::InitHistograms()
+{
+ gStyle->SetOptStat();
+
+ fHistResidualX.clear();
+ fHistResidualY.clear();
+ fHistResidualT.clear();
+ fHistPullX.clear();
+ fHistPullY.clear();
+ fHistPullT.clear();
-// --------------------------------------------------------------------------------------------------------------------
+ fHistPointsPerHit.clear();
+ fHistHitsPerPoint.clear();
+
+ fHistEfficiencyR.clear();
+ fHistEfficiencyXY.clear();
+
+ const int nStations = fpDetectorInterface->GetNtrackingStations();
+
+ fHistResidualX.reserve(nStations);
+ fHistResidualY.reserve(nStations);
+ fHistResidualT.reserve(nStations);
+
+ fHistPullX.reserve(nStations);
+ fHistPullY.reserve(nStations);
+ fHistPullT.reserve(nStations);
+
+ fHistPointsPerHit.reserve(nStations);
+ fHistHitsPerPoint.reserve(nStations);
+
+ fHistEfficiencyR.reserve(nStations);
+ fHistEfficiencyXY.reserve(nStations);
+
+ // Create histograms
+ TString histName = "";
+ TString histTitle = "";
+ int nBins = -1;
+ double rMin = 0.;
+ double rMax = 0.;
+ for (int iSt = 0; iSt < nStations; ++iSt) {
+ // Residuals
+ histName = Form("h1_sts_ResidualX_st%d", iSt);
+ histTitle = Form("STS: Station %d: Residual X; (X_{reco} - X_{MC}) [cm]", iSt);
+ std::tie(nBins, rMin, rMax) = fRangeResidualX;
+ fHistResidualX.emplace_back(histName, histTitle, nBins, rMin, rMax);
+
+ histName = Form("h1_sts_ResidualY_st%d", iSt);
+ histTitle = Form("STS: Station %d: Residual Y; (Y_{reco} - Y_{MC}) [cm]", iSt);
+ std::tie(nBins, rMin, rMax) = fRangeResidualY;
+ fHistResidualY.emplace_back(histName, histTitle, nBins, rMin, rMax);
+
+ histName = Form("h1_sts_ResidualT_st%d", iSt);
+ histTitle = Form("STS: Station %d: Residual T; (T_{reco} - T_{MC}) [ns]", iSt);
+ std::tie(nBins, rMin, rMax) = fRangeResidualT;
+ fHistResidualT.emplace_back(histName, histTitle, nBins, rMin, rMax);
+
+ // Pulls
+ histName = Form("h1_sts_PullX_st%d", iSt);
+ histTitle = Form("STS: Station %d: Pull X; (X_{reco} - X_{MC}) / #sigmaX_{reco}", iSt);
+ std::tie(nBins, rMin, rMax) = fRangePullX;
+ fHistPullX.emplace_back(histName, histTitle, nBins, rMin, rMax);
+
+ histName = Form("h1_sts_PullY_st%d", iSt);
+ histTitle = Form("STS: Station %d: Pull Y; (Y_{reco} - Y_{MC}) / #sigmaY_{reco}", iSt);
+ std::tie(nBins, rMin, rMax) = fRangePullY;
+ fHistPullY.emplace_back(histName, histTitle, nBins, rMin, rMax);
+
+ histName = Form("h1_sts_PullT_st%d", iSt);
+ histTitle = Form("STS: Station %d: Pull T; (T_{reco} - T_{MC}) / #sigmaT_{reco}", iSt);
+ std::tie(nBins, rMin, rMax) = fRangePullT;
+ fHistPullT.emplace_back(histName, histTitle, nBins, rMin, rMax);
+
+ // MC points per one hit
+ histName = Form("h1_sts_PointsPerHit_st%d", iSt);
+ histTitle = Form("STS: Station %d: MC Points per Hit; N_{MC Points} per hit", iSt);
+ std::tie(nBins, rMin, rMax) = fRangePointsPerHit;
+ fHistPointsPerHit.emplace_back(histName, histTitle, nBins, rMin, rMax);
+
+ // Hits per one mc point
+ histName = Form("h1_sts_HitsPerPoint_st%d", iSt);
+ histTitle = Form("STS: Station %d: Hit per MC points; N_{hits} per MC point", iSt);
+ std::tie(nBins, rMin, rMax) = fRangeHitsPerPoint;
+ fHistHitsPerPoint.emplace_back(histName, histTitle, nBins, rMin, rMax);
+
+ // Efficiencies
+ double xMax = fpDetectorInterface->GetXmax(iSt); /// TODO: test ranges
+ double yMax = fpDetectorInterface->GetYmax(iSt); /// TODO: test ranges
+
+ // Efficiency vs. distance of point from center
+ histName = Form("pr1_sts_EfficiencyR_st%d", iSt);
+ histTitle = Form("STS: Station %d: Efficiency R; R [cm]", iSt);
+ rMin = 0.;
+ rMax = sqrt(xMax * xMax + yMax * yMax);
+ nBins = 100;
+ fHistEfficiencyR.emplace_back(histName, histTitle, nBins, rMin, rMax);
+ fHistEfficiencyR[iSt].SetOptStat(1110);
+
+ // Efficiency vs. x and y
+ histName = Form("pr1_sts_EfficiencyXY_st%d", iSt);
+ histTitle = Form("STS: Station %d: Efficiency XY; X [cm]; Y [cm]", iSt);
+ nBins = 50;
+ fHistEfficiencyXY.emplace_back(histName, histTitle, nBins, -xMax, xMax, nBins, -yMax, yMax);
+ fHistEfficiencyXY[iSt].SetOptStat(10);
+ }
+
+ // Register histograms
+ for (int iSt = 0; iSt < nStations; ++iSt) {
+ RegisterHist(&fHistResidualX[iSt]);
+ RegisterHist(&fHistResidualY[iSt]);
+ RegisterHist(&fHistResidualT[iSt]);
+
+ RegisterHist(&fHistPullX[iSt]);
+ RegisterHist(&fHistPullY[iSt]);
+ RegisterHist(&fHistPullT[iSt]);
+
+ RegisterHist(&fHistPointsPerHit[iSt]);
+ RegisterHist(&fHistHitsPerPoint[iSt]);
+
+ RegisterHist(&fHistEfficiencyXY[iSt]);
+ RegisterHist(&fHistEfficiencyR[iSt]);
+ }
+
+ return kSUCCESS;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+CbmMatch CbmTrackingInputQaSts::MatchHits(const CbmStsHit* pHit, int iHit)
+{
+ CbmMatch res; // Matching result
+
+ // Front and back cluster indeces
+ const int iClusterF = pHit->GetFrontClusterId();
+ if (iClusterF < 0) {
+ LOG(error) << "STD: hit (id = " << iHit << ") has incorrect front cluster index: " << iClusterF;
+ throw std::runtime_error("STS tracking input QA: wrong front cluster index");
+ }
+
+ const int iClusterB = pHit->GetBackClusterId();
+ if (iClusterB < 0) {
+ LOG(error) << "STD: hit (id = " << iHit << ") has incorrect back cluster index: " << iClusterF;
+ throw std::runtime_error("STS tracking input QA: wrong back cluster index");
+ }
+
+ // Front and back clusters of a hit
+ const auto* pClusterF = dynamic_cast<const CbmStsCluster*>(fpClusters->At(iClusterF));
+ if (!pClusterF) {
+ LOG(error) << "STD: front cluster does not exist for hit (id = " << iHit << ')';
+ throw std::runtime_error("STS tracking input QA: front cluster not found");
+ }
+
+ const auto* pClusterB = dynamic_cast<const CbmStsCluster*>(fpClusters->At(iClusterB));
+ if (!pClusterB) {
+ LOG(error) << "STD: back cluster does not exist for hit (id = " << iHit << ')';
+ throw std::runtime_error("STS tracking input QA: back cluster not found");
+ }
+
+ const CbmMatch* pClusterMatchF = dynamic_cast<const CbmMatch*>(fpClusterMatches->At(iClusterF));
+ const CbmMatch* pClusterMatchB = dynamic_cast<const CbmMatch*>(fpClusterMatches->At(iClusterB));
+
+ // Check addresses
+ const int addressHit = pHit->GetAddress();
+ const int addressClusterF = pClusterF->GetAddress();
+ const int addressClusterB = pClusterB->GetAddress();
+ if (addressHit != addressClusterF || addressHit != addressClusterB) {
+ LOG(error) << "STS hit (id = " << iHit << ") and its front and(or) back clusters has different addresses: "
+ << "hit address = " << addressHit << ", front cluster address = " << addressClusterF
+ << ", back cluster address = " << addressClusterB;
+ throw std::runtime_error("STS tracking input QA: inconsistent in hit and clusters addresses");
+ }
+
+ // CUSTOM MATCHING: we choose only those links, which are presented both in front and back clusters
+ // The same matching is used in L1 tracking
+ 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) {
+ res.AddLink(linkF);
+ res.AddLink(linkB);
+ }
+ }
+ }
+
+ return res; // Rely on NRVO
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+double CbmTrackingInputQaSts::ParticleMass(int pdg)
+{
+ 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";
+
+ return 0.;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+InitStatus CbmTrackingInputQaSts::ReadAndCreateDataBranches()
+{
+ auto* pManager = FairRootManager::Instance();
+ if (!pManager) {
+ LOG(fatal) << "FairRootManager was not found";
+ }
+
+ fpDigiManager = CbmDigiManager::Instance();
+ fpDigiManager->Init(); // NOTE: is initialized only once
+
+ // ************************************************
+ // ** Reconstructed data branches initialization **
+ // ************************************************
+
+ fpTimeSlice = dynamic_cast<CbmTimeSlice*>(pManager->GetObject("TimeSlice."));
+ if (!fpTimeSlice) {
+ LOG(error) << "Time slice was not found";
+ return kERROR;
+ }
+
+ fpHits = dynamic_cast<TClonesArray*>(FairRootManager::Instance()->GetObject("StsHit"));;
+ if (!fpHits) {
+ LOG(error) << "Hits input array was not found for STS";
+ return kERROR;
+ }
+
+ fpClusters = dynamic_cast<TClonesArray*>(FairRootManager::Instance()->GetObject("StsCluster"));
+ if (!fpClusters) {
+ LOG(error) << "Cluster input array was not found for STS";
+ return kERROR;
+ }
+
+ // MC data initialization
+ fpMcManager = dynamic_cast<CbmMCDataManager*>(pManager->GetObject("MCDataManager"));
+ fIsMcPresent = bool(fpMcManager);
+
+ if (fIsMcPresent) {
+ fpMcEventList = dynamic_cast<CbmMCEventList*>(pManager->GetObject("MCEventList."));
+ fpMcTracks = fpMcManager->InitBranch("MCTrack");
+ fpMcPoints = fpMcManager->InitBranch("StsPoint");
+ fpHitMatches = dynamic_cast<TClonesArray*>(FairRootManager::Instance()->GetObject("StsHitMatch"));
+ fpClusterMatches = dynamic_cast<TClonesArray*>(FairRootManager::Instance()->GetObject("StsClusterMatch"));
+
+ if (!fpMcEventList) {
+ LOG(error) << "MCEventList data were not found";
+ return kERROR;
+ }
+
+ if (!fpMcTracks) {
+ LOG(error) << "MC tracks data were not found";
+ return kERROR;
+ }
+
+ if (!fpMcPoints) {
+ LOG(error) << "MC points data were not found for STS";
+ return kERROR;
+ }
+
+ if (!fpHitMatches) {
+ LOG(error) << "Hits match data were not found for STS";
+ return kERROR;
+ }
+ }
+
+ return kSUCCESS;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void CbmTrackingInputQaSts::ResolutionQa()
+{
+ // Resolution QA is impossible without MC information
+ if (!fIsMcPresent) { return; }
+
+ const int nHits = fpHits->GetEntriesFast(); // Number of hits in event/ts
+ const int nStations = fpDetectorInterface->GetNtrackingStations(); // Number of stations in event/ts
+ const int nMcEvents = fpMcEventList->GetNofEvents(); // Number of MC events in event/ts
+
+ std::vector<std::vector<int>> nHitsPerMcPoint; // Number of hits per one MC point in different MC events
+ nHitsPerMcPoint.resize(nMcEvents);
+
+ // ** Loop over MC events within event/ts **
+ for (int iE = 0; iE < nMcEvents; ++iE) {
+ const int fileId = fpMcEventList->GetFileIdByIndex(iE);
+ const int eventId = fpMcEventList->GetEventIdByIndex(iE);
+ const int nPoints = fpMcPoints->Size(fileId, eventId);
+ nHitsPerMcPoint[iE].resize(nPoints, 0);
+ } // Loop over MC events within event/ts: end
+
+
+ // ** Loop over hits within event/ts **
+ for (int iHit = 0; iHit < nHits; ++iHit) {
+ const auto* pHit = dynamic_cast<const CbmStsHit*>(fpHits->At(iHit));
+ if (!pHit) {
+ LOG(error) << "STS hit with index " << iHit << " does not exist in \"StsHit\" array";
+ throw std::runtime_error("STS tracking input QA: hit does not exist");
+ }
+
+ // Check station index of a hit
+ const int iStation = fpDetectorInterface->GetTrackingStationIndex(pHit);
+ if (iStation < 0 || iStation >= nStations) {
+ LOG(error) << "STS hit with index " << iHit << ": tracking station index = "
+ << iStation << " is out of range [0, " << nStations << ']';
+ throw std::runtime_error("STS tracking input QA: hit has inconsistent station index");
+ }
+
+ // Get custom match for the hit
+ const CbmMatch hitMatch = MatchHits(pHit, iHit); // throws std::runtime_error
+
+ // Fill number of hits per one MC point vector and update the number of hits per a given MC point
+ int nMcPoints = 0; // Number of non-noisy MC points per event
+ for (int iLink = 0; iLink < hitMatch.GetNofLinks(); ++iLink) {
+ const auto& link = hitMatch.GetLink(iLink);
+ if (link.GetIndex() >= 0) { // Select only non-noisy links (non-noisy digis)
+ ++nMcPoints;
+ const int iE = fpMcEventList->GetEventIndex(link);
+ if (iE < 0 || iE >= nMcEvents) {
+ LOG(error) << "STS: link points to a non-existing MC event (iHit = " << iHit << ')';
+ throw std::runtime_error("STS tracking input QA: link points to non-existing MC event");
+ }
+ if (link.GetIndex() >= int(nHitsPerMcPoint[iE].size())) {
+ LOG(error) << "STS: link points to a non-existing MC point (iHit = " << iHit << ')';
+ throw std::runtime_error("STS tracking input QA: link points to non-existing MC point");
+ }
+ nHitsPerMcPoint[iE][link.GetIndex()]++;
+ }
+ }
+
+ fHistPointsPerHit[iStation].Fill(nMcPoints);
+
+ // Select only hits with one MC point
+ if (nMcPoints != 1) { continue; }
+
+ // Take the best link in the match
+ const auto& bestLink = hitMatch.GetMatchedLink();
+
+ // Skip noise
+ if (bestLink.GetIndex() < 0) { continue; }
+
+ // Get MC point
+ const auto* pPoint = dynamic_cast<const CbmStsPoint*>(fpMcPoints->Get(bestLink));
+ if (!pPoint) {
+ LOG(error) << "STS: MC point does not exist (iHit = " << iHit << ')';
+ throw std::runtime_error("STS tracking input QA: MC point does not exist"); // NOTE:
+ }
+
+ // TODO: Study the distributions of |pointZIn - stationZ|, |pointZout - stationZ| and |hitZ - stationZ|
+ //const double stationZ = fpDetectorInterface->GetZ(iStation);
+
+ // ***********************************
+ // ** Calculate residuals and pools **
+ // ***********************************
+
+ // Get MC event time
+ double t0 = fpMcEventList->GetEventTime(bestLink);
+ if (t0 < 0) {
+ LOG(error) << "STS: MC event time is undefined (iHit = " << iHit << ')';
+ throw std::runtime_error("STS tracking input QA: the MC event time is undefined");
+ }
+
+ // Get time, space position and momenta components for the MC point (values are taken for the "In" point)
+ double mcX = pPoint->GetXIn(); // [cm]
+ double mcY = pPoint->GetYIn(); // [cm]
+ double mcZ = pPoint->GetZIn(); // [cm]
+ double mcT = pPoint->GetTime() + t0; // [ns]
+ double mcPx = pPoint->GetPx(); // [GeV/c]
+ double mcPy = pPoint->GetPy(); // [GeV/c]
+ double mcPz = pPoint->GetPz(); // [GeV/c]
+
+ // Skip slow particles (TODO: why pZ, not p?
+ if (fabs(pPoint->GetPzOut()) < fMinMomentum) { continue; }
+
+ // Difference between z components of MC in point and the hit
+ double dz = pHit->GetZ() - mcZ; // [cm]
+
+ // Propagate MC-point x and y "In" coordinates to z-plane of the hit
+ mcX += dz * mcPx / mcPz;
+ mcY += dz * mcPy / mcPz;
+
+ // Propagete MC-point "In" time to z-plane of the hit
+ int pdgCode = pPoint->GetPid();
+ double mass = ParticleMass(pdgCode);
+ constexpr double speedOfLight {29.9792458}; // cm/ns
+ TVector3 mom;
+
+ mcT += dz / (mcPz * speedOfLight) * sqrt(mass * mass + mom.Mag2());
+
+ double dx = pHit->GetX() - mcX;
+ double dy = pHit->GetY() - mcY;
+ double dt = pHit->GetTime() - mcT;
+ double rmsX = pHit->GetDx();
+ double rmsY = pHit->GetDy();
+ double rmsT = pHit->GetTimeError();
+
+ fHistResidualX[iStation].Fill(dx);
+ fHistResidualY[iStation].Fill(dy);
+ fHistResidualT[iStation].Fill(dt);
+
+ fHistPullX[iStation].Fill(dx / rmsX);
+ fHistPullY[iStation].Fill(dy / rmsY);
+ fHistPullT[iStation].Fill(dt / rmsT);
+ }// Loop over hits within event/ts: end
+
+ // *********************
+ // ** Fill efficiency **
+ // *********************
+
+ // ** Loop over MC events within event/ts **
+ for (int iE = 0; iE < nMcEvents; ++iE) {
+ const int fileId = fpMcEventList->GetFileIdByIndex(iE);
+ const int eventId = fpMcEventList->GetEventIdByIndex(iE);
+ const int nPoints = fpMcPoints->Size(fileId, eventId);
+
+ // ** Loop over MC points **
+ for (int iP = 0; iP < nPoints; ++iP) {
+ const auto* pPoint = dynamic_cast<CbmStsPoint*>(fpMcPoints->Get(fileId, eventId, iP));
+ if (!pPoint) {
+ LOG(error) << "MC point does not exist for iE = " << iE << ", iP = " << iP;
+ throw std::runtime_error("STS tracking input QA: MC point does not exist");
+ }
+
+ int address = pPoint->GetDetectorID();
+ int iStation = fpDetectorInterface->GetTrackingStationIndex(address);
+ if (iStation < 0 || iStation >= nStations) {
+ LOG(error) << "STS MC point with index " << iP << ": tracking station index = "
+ << iStation << " is out of range [0, " << nStations << ']';
+ throw std::runtime_error("STS tracking input QA: MC point has inconsistent station index");
+ }
+ fHistHitsPerPoint[iStation].Fill(nHitsPerMcPoint[iE][iP]);
+
+ double pointDistance = sqrt(pPoint->GetXIn() * pPoint->GetXIn() + pPoint->GetYIn() * pPoint->GetYIn()); // [cm]
+ fHistEfficiencyR[iStation].Fill(pointDistance, (nHitsPerMcPoint[iE][iP] > 0));
+ fHistEfficiencyXY[iStation].Fill(pPoint->GetXIn(), pPoint->GetYIn(), (nHitsPerMcPoint[iE][iP] > 0));
+ }
+
+ } // Loop over MC events within event/ts: end
+}
+
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void CbmTrackingInputQaSts::SetParContainers()
+{
+ std::cout << "\033[1;32mCALL:\033[0mCbmTrackingInputQaSts::SetParContainers()\n";
+}
-// --------------------------------------------------------------------------------------------------------------------
diff --git a/reco/L1/qa/CbmTrackingInputQaSts.h b/reco/L1/qa/CbmTrackingInputQaSts.h
index f866655a99..7a0b324baa 100644
--- a/reco/L1/qa/CbmTrackingInputQaSts.h
+++ b/reco/L1/qa/CbmTrackingInputQaSts.h
@@ -4,52 +4,242 @@
/***************************************************************************************************
* @file CbmTrackingInputQaSts.h
- * @brief Class for the tracking input QA for STS (declaration)
- * @since 27.06.2022
+ * @brief Base class for the tracking input QA (declaration)
+ * @since 24.06.2022
* @author S.Zharko <s.zharko@gsi.de>
***************************************************************************************************/
#ifndef CbmTrackingInputQaSts_h
#define CbmTrackingInputQaSts_h 1
-#include "CbmTrackingInputQaBase.h"
+#include "CbmMCDataManager.h"
+#include "CbmQaCanvas.h"
+#include "CbmQaHist.h"
+#include "CbmSetup.h"
+#include "CbmTrackingDetectorInterfaceBase.h" // Communicate via tracking detector interface
-#include "FairRootManager.h"
+#include <FairTask.h>
-#include "TClonesArray.h"
+#include "TFolder.h"
+#include "TH1D.h"
+#include "TH1I.h"
+#include "TH2D.h"
+#include "TParameter.h"
+#include "TProfile.h"
+#include "TString.h"
+#include <tuple>
+#include <memory>
+#include <vector>
+
+class CbmDigiManager;
+class CbmMCEventList;
class CbmMCDataArray;
-class CbmMatch;
+class CbmTimeSlice;
+class TClonesArray;
class CbmCluster;
+class CbmMatch;
+class CbmStsHit;
-class CbmTrackingInputQaSts : public CbmTrackingInputQaBase {
+/// Class CbmTrackingInputQaSts is an abstract FairTask inherited class, providing interface for a
+/// tracking detector QA including input geometry and data checks. One should inherit this class
+/// for a particular detector subsystem.
+///
+class CbmTrackingInputQaSts : public FairTask {
public:
- CbmTrackingInputQaSts(int verbose = 1);
+ /// Constructor from the task name and verbosity level
+ CbmTrackingInputQaSts(int verbosity = 0);
+
+ /// Destructor
~CbmTrackingInputQaSts();
+ /// Suppress copy and move
+ CbmTrackingInputQaSts(const CbmTrackingInputQaSts&) = delete;
+ CbmTrackingInputQaSts(CbmTrackingInputQaSts&&) = delete;
+ CbmTrackingInputQaSts& operator=(const CbmTrackingInputQaSts&) = delete;
+ CbmTrackingInputQaSts& operator=(CbmTrackingInputQaSts&&) = delete;
+
+ /// TTask: Executes the task (processes a timeslice)
+ /// \param option Specify action for the TTask::Exec (not defined yet)
+ void Exec(Option_t* /*option*/);
+
+ /// FairTask: Action at hte end of the run
+ void Finish();
+
+ /// Gets maximum allowed distance between z-components of hit/MC-point position and the station center [cm]
+ double GetMaxDistanceZ() const { return fMaxDistanceZ; }
+
+ /// Gets maximum allowed difference between pulls distributions RMS and the expected RMS value of RMS = 1
+ double GetMaxPullsRmsDiff() const { return fMaxPullsRmsDiff; }
+
+ /// Gets minimum particle momentum
+ double GetMinMomentum() const { return fMinMomentum; }
+
+ /// Fits and draws histograms
+ TFolder& GetQa();
+
+ /// FairTask: Task initialization an the beginning of the run
InitStatus Init();
+ /// FairTask: Task reinitialization
+ InitStatus ReInit() { return Init(); }
- /// Gets pointer to hits array of a particular detector type
- TClonesArray* GetHitsInput() { return dynamic_cast<TClonesArray*>(FairRootManager::Instance()->GetObject("StsHit")); }
+ /// Sets maximum allowed distance between z-components of hit/MC-point position and the station center [cm]
+ void SetMaxDistanceZ(double dist) { fMaxDistanceZ = dist; }
- /// Gets pointer to clusters array of a particular detector type
- TClonesArray* GetClustersInput()
- {
- return dynamic_cast<TClonesArray*>(FairRootManager::Instance()->GetObject("StsCluster"));
- }
+ /// Sets maximum allowed difference between pulls distributions RMS and the expected RMS value of RMS = 1
+ void SetMaxPullsRmsDiff(double rmsMaxDiff) { fMaxPullsRmsDiff = rmsMaxDiff; }
+
+ /// Sets minimum particle momentum
+ void SetMinMomentum(double pMin) { fMinMomentum = pMin; }
+
+ /// FairTask: Set parameter containers (if any)
+ void SetParContainers();
+
+ ClassDef(CbmTrackingInputQaSts, 0); // Base class for tracking input QA
+
+private:
+ // ***************
+ // ** Functions **
+ // ***************
+
+ /// Checks accumulated distributions
+ bool CheckDistributions();
+
+ /// Deinitialize data
+ void DeInit();
+
+ /// Fills distributions of residuals
+ bool FillHistosResiduals();
- /// Gets pointer to hit matches array for a particular detector type
- TClonesArray* GetHitMatchesInput()
+ /// Checks geometry
+ InitStatus GeometryQa();
+
+ /// Returns particle mass [GeV/c2] by its pdg code
+ static double ParticleMass(int pdg);
+
+ /// Initializes canvases
+ InitStatus InitCanvases();
+
+ /// Initializes base histograms
+ InitStatus InitHistograms();
+
+ /// Match hits with MC points and fills out the vector
+ /// \param pHit A pointer to the hit to be matched
+ /// \param iHit An index of the hit in the hits array
+ /// \return Match object
+ CbmMatch MatchHits(const CbmStsHit* pHit, int iHit);
+
+ /// Initializes input data branches
+ InitStatus ReadAndCreateDataBranches();
+
+ /// Registers histogram
+ /// \param pHisto
+ void RegisterHist(TH1* pHisto)
{
- return dynamic_cast<TClonesArray*>(FairRootManager::Instance()->GetObject("StsHitMatch"));
+ fHistList.push_back(pHisto);
+ pHisto->SetDirectory(0);
}
- /// Gets pointer to MC points for a particular detector type
- CbmMCDataArray* GetMcPointsInput() { return fpMcManager->InitBranch("StsPoint"); }
+ /// Checks resolution
+ /// A general scheme of resolution checks is provided in this method. One should specify particular operations in the
+ /// derived classes for each particular detector....
+ void ResolutionQa();
+
+ // *****************
+ // ** Data fields **
+ // *****************
+
+ // ** Flags and cuts **
+ bool fIsMcPresent {false}; ///< Flag: true - all procedures involving MC will be processed, otherwice data will be processed only
+
+ bool fIsQaPassed {true}; ///< Flag: true - QA is successfull, false - Qa is failed
+
+ /// Maximum allowed distance between z-components of hit/MC-point position and the station center [cm]
+ double fMaxDistanceZ {1.0};
+
+ double fMinMomentum {0.01}; ///< Minimum value of particle momentum z component [GeV/c]
+
+ double fMaxPullsRmsDiff {0.05}; ///< Maximum difference of the pulls RMS from unity
+
+ // ** Input arrays **
+
+ CbmTrackingDetectorInterfaceBase* fpDetectorInterface {nullptr}; ///< Pointer to current tracking detector I/F
+ CbmTimeSlice* fpTimeSlice {nullptr}; ///< Pointer to current time slice
+ CbmDigiManager* fpDigiManager {nullptr}; ///< Pointer to digi manager
+
+ CbmMCEventList* fpMcEventList {nullptr};
+ CbmMCDataManager* fpMcManager {nullptr};
+
+ TClonesArray* fpClusters {nullptr}; ///< Clusters
+ TClonesArray* fpHits {nullptr}; ///< Hits
+ TClonesArray* fpClusterMatches {nullptr}; ///< Matches for clusters (used only in STS)
+ TClonesArray* fpHitMatches {nullptr}; ///< Matches for hits, provided with CbmMatchRecoToMC task
+
+ CbmMCDataArray* fpMcTracks {nullptr};
+ CbmMCDataArray* fpMcPoints {nullptr};
+
+
+ // ** Output folders **
+
+ TFolder fOutFolder {"TrackingInputQaSts", "TrackingInputQaSts"}; ///< Output folder for this QA task
+ TFolder* fpOutFolderHists {nullptr}; ///< Subfolder for raw histograms
+
+ std::vector<TH1*> fHistList; ///< List of the pointers to all histograms contained in the class
+
+ // ** Histograms **
+
+ TParameter<long> fNevents {"nEvents", 0}; ///< Number of processed events
+
+ std::vector<CbmQaHist<TH1D>> fHistResidualX; ///< Residual distributions for X vs index of a station
+ std::vector<CbmQaHist<TH1D>> fHistResidualY; ///< Residual distributions for X vs index of a station
+ std::vector<CbmQaHist<TH1D>> fHistResidualT; ///< Residual distributions for T vs index of a station
+
+ std::vector<CbmQaHist<TH1D>> fHistPullX; ///< Pulls distribution for X vs. index of a station
+ std::vector<CbmQaHist<TH1D>> fHistPullY; ///< Pulls distribution for Y vs. index of a station
+ std::vector<CbmQaHist<TH1D>> fHistPullT; ///< Pulls distribution for T vs. index of a station
+
+ std::vector<CbmQaHist<TH1D>> fHistPointsPerHit; ///< Distribution of MC points per one hit vs. index of a station
+ std::vector<CbmQaHist<TH1D>> fHistHitsPerPoint; ///< Distribution of hits per one MC point vs. index of a station
+
+ // TODO: There is a class TEfficiency in the ROOT, probably it's a good idea to use it here (S.Zharko)
+ std::vector<CbmQaHist<TProfile2D>> fHistEfficiencyXY; ///<
+ std::vector<CbmQaHist<TProfile>> fHistEfficiencyR; ///<
+
+ static constexpr std::tuple<int, double, double> fRangeResidualX {100, -0.01, 0.01};
+ static constexpr std::tuple<int, double, double> fRangeResidualY {100, -0.03, 0.03};
+ static constexpr std::tuple<int, double, double> fRangeResidualT {100, -25., 25.};
+
+ static constexpr std::tuple<int, double, double> fRangePullX {100, -5., 5.};
+ static constexpr std::tuple<int, double, double> fRangePullY {100, -5., 5.};
+ static constexpr std::tuple<int, double, double> fRangePullT {100, -5., 5.};
+
+ static constexpr std::tuple<int, double, double> fRangePointsPerHit {10, -0.5, 9.5};
+ static constexpr std::tuple<int, double, double> fRangeHitsPerPoint {10, -0.5, 9.5};
+
+
+ // ** Canvases **
+
+ CbmQaCanvas fCanvResidualX {"c_sts_residualX", "STS X Residual distributions", 1500, 1000 };
+ CbmQaCanvas fCanvResidualY {"c_sts_residualY", "STS Y Residual distributions", 1500, 1000 };
+ CbmQaCanvas fCanvResidualT {"c_sts_residualT", "STS T Residual distributions", 1500, 1000 };
+
+ CbmQaCanvas fCanvPullX {"c_sts_pullX", "STS X Pull distributions", 1500, 1000};
+ CbmQaCanvas fCanvPullY {"c_sts_pullY", "STS Y Pull distributions", 1500, 1000};
+ CbmQaCanvas fCanvPullT {"c_sts_pullT", "STS T Pull distributions", 1500, 1000};
+
+ CbmQaCanvas fCanvPointsPerHit {"c_sts_PointsPerHit", "STS MC Points per Hit", 1500, 1000};
+ CbmQaCanvas fCanvHitsPerPoint {"c_sts_HitsPerPoint", "STS Hits per MC Point", 1500, 1000};
+
+ CbmQaCanvas fCanvEfficiencyR {"c_sts_EfficiencyR", "STS Efficiency vs. R", 1500, 1000};
+ CbmQaCanvas fCanvEfficiencyXY {"c_sts_EfficiencyXY", "STS Efficiency vs. X and Y", 1500, 1000};
+
+
+ // ************
+ // ** Output **
+ // ************
- ClassDef(CbmTrackingInputQaSts, 0);
};
-#endif // CbmTrackingInputQaSts_h
+#endif
--
GitLab