From 3767a70832e38d8f8c6fc2300f2becb37e67507c Mon Sep 17 00:00:00 2001
From: "s.zharko@gsi.de" <s.zharko@gsi.de>
Date: Wed, 28 Jun 2023 20:30:41 +0200
Subject: [PATCH] CA: Config and uniform distribution for ideal hit finder
---
core/base/CbmMatchRecoToMC.cxx | 30 +--
core/base/CbmMatchRecoToMC.h | 8 +-
core/base/CbmTrackingDetectorInterfaceBase.h | 2 +-
core/data/tof/CbmTofHit.cxx | 1 +
macro/CMakeLists.txt | 1 +
macro/L1/CMakeLists.txt | 1 +
macro/L1/configs/CMakeLists.txt | 4 +
macro/L1/configs/README | 8 +
macro/L1/configs/config_ideal_hits_mcbm.yaml | 42 ++++
macro/mcbm/mcbm_reco_event.C | 17 +-
macro/qa/qa_compare_ca.C | 45 ++++
macro/run/run_reco.C | 34 ---
reco/L1/L1Algo/utils/CaAlgoRandom.h | 20 +-
reco/L1/utils/CbmCaIdealHitProducer.cxx | 42 +---
reco/L1/utils/CbmCaIdealHitProducer.h | 21 +-
reco/L1/utils/CbmCaIdealHitProducerDet.h | 231 ++++++++++++-------
16 files changed, 314 insertions(+), 193 deletions(-)
create mode 100644 macro/L1/CMakeLists.txt
create mode 100644 macro/L1/configs/CMakeLists.txt
create mode 100644 macro/L1/configs/README
create mode 100644 macro/L1/configs/config_ideal_hits_mcbm.yaml
create mode 100644 macro/qa/qa_compare_ca.C
diff --git a/core/base/CbmMatchRecoToMC.cxx b/core/base/CbmMatchRecoToMC.cxx
index a01144afa2..3e11ab73f4 100644
--- a/core/base/CbmMatchRecoToMC.cxx
+++ b/core/base/CbmMatchRecoToMC.cxx
@@ -118,7 +118,7 @@ InitStatus CbmMatchRecoToMC::Init()
ReadAndCreateDataBranches();
// Notification to a user about matching suppression
- if (fbSuppressMatching) {
+ if (fbSuppressReMatching) {
std::stringstream msg;
msg << "\033[1;31mCbmMatchRecoToMC (\"" << fName << "\"): the cluster and hit matching routines for MVD, STS, ";
msg << "MuCh, TRD, TOF will be suppressed by a request from CbmMatchRecoToMC::SuppressHitMatching():\033[0m\n";
@@ -130,7 +130,7 @@ InitStatus CbmMatchRecoToMC::Init()
void CbmMatchRecoToMC::Exec(Option_t* /*opt*/)
{
- if (!fbSuppressMatching) {
+ if (!fbSuppressReMatching) {
if (fMvdHitMatches != nullptr) fMvdHitMatches->Delete();
if (fMvdClusterMatches != nullptr) fMvdClusterMatches->Delete();
if (fStsClusterMatches != nullptr) fStsClusterMatches->Delete();
@@ -149,7 +149,7 @@ void CbmMatchRecoToMC::Exec(Option_t* /*opt*/)
// ----- MVD -----
- if (!fbSuppressMatching) {
+ if (!fbSuppressReMatching) {
// MVD: digi->hit
if (fMvdHits && fMvdHitMatches && !fMvdCluster) { MatchHitsMvd(fMvdHits, fMvdHitMatches); }
else {
@@ -161,7 +161,7 @@ void CbmMatchRecoToMC::Exec(Option_t* /*opt*/)
}
// ----- STS -----
- if (!fbSuppressMatching) {
+ if (!fbSuppressReMatching) {
// STS: digi->cluster
MatchClusters(ECbmModuleId::kSts, fStsClusters, fStsClusterMatches);
// STS: cluster->hit
@@ -171,7 +171,7 @@ void CbmMatchRecoToMC::Exec(Option_t* /*opt*/)
MatchStsTracks(fMvdHitMatches, fStsHitMatches, fMvdPoints, fStsPoints, fStsTracks, fStsTrackMatches);
// ----- MUCH -----
- if (!fbSuppressMatching) {
+ if (!fbSuppressReMatching) {
// MUCH: digi->cluster
MatchClusters(ECbmModuleId::kMuch, fMuchClusters, fMuchClusterMatches);
// MUCH: cluster->hit
@@ -188,19 +188,19 @@ void CbmMatchRecoToMC::Exec(Option_t* /*opt*/)
// TRD
if (fTrdClusters && fTrdHits) { // MC->digi->cluster->hit->track
- if (!fbSuppressMatching) {
+ if (!fbSuppressReMatching) {
MatchClusters(ECbmModuleId::kTrd, fTrdClusters, fTrdClusterMatches);
MatchHits(fTrdClusterMatches, fTrdHits, fTrdHitMatches);
}
MatchTracks(fTrdHitMatches, fTrdPoints, fTrdTracks, fTrdTrackMatches);
}
else if (fTrdHits) { // MC->hit->track
- if (!fbSuppressMatching) { MatchHitsToPoints(fTrdPoints, fTrdHits, fTrdHitMatches); }
+ if (!fbSuppressReMatching) { MatchHitsToPoints(fTrdPoints, fTrdHits, fTrdHitMatches); }
MatchTracks(fTrdHitMatches, fTrdPoints, fTrdTracks, fTrdTrackMatches);
}
// TOF: (Digi->MC)+(Hit->Digi)=>(Hit->MC)
- if (!fbSuppressMatching) { MatchHitsTof(fTofHitDigiMatches, fTofHits, fTofHitMatches); }
+ if (!fbSuppressReMatching) { MatchHitsTof(fTofHitDigiMatches, fTofHits, fTofHitMatches); }
//static Int_t eventNo = 0;
LOG(info) << "CbmMatchRecoToMC::Exec eventNo=" << fEventNumber++;
}
@@ -236,7 +236,7 @@ void CbmMatchRecoToMC::ReadAndCreateDataBranches()
if (nullptr != fStsClusters) {
fStsClusterMatches = static_cast<TClonesArray*>(ioman->GetObject("StsClusterMatch"));
if (nullptr == fStsClusterMatches) {
- fbSuppressMatching = false; // Never the less, do matching, because there are no matches
+ fbSuppressReMatching = false; // Ensure, that matching is done, if matches are absent
fStsClusterMatches = new TClonesArray("CbmMatch", 100);
ioman->Register("StsClusterMatch", "STS", fStsClusterMatches, IsOutputBranchPersistent("StsClusterMatch"));
}
@@ -246,7 +246,7 @@ void CbmMatchRecoToMC::ReadAndCreateDataBranches()
if (nullptr != fStsHits) {
fStsHitMatches = static_cast<TClonesArray*>(ioman->GetObject("StsHitMatch"));
if (nullptr == fStsHitMatches) {
- fbSuppressMatching = false; // Never the less, do matching, because there are no matches
+ fbSuppressReMatching = false; // Ensure, that matching is done, if matches are absent
fStsHitMatches = new TClonesArray("CbmMatch", 100);
ioman->Register("StsHitMatch", "STS", fStsHitMatches, IsOutputBranchPersistent("StsHitMatch"));
}
@@ -280,7 +280,7 @@ void CbmMatchRecoToMC::ReadAndCreateDataBranches()
if (nullptr != fTrdClusters) {
fTrdClusterMatches = static_cast<TClonesArray*>(ioman->GetObject("TrdClusterMatch"));
if (nullptr == fTrdClusterMatches) {
- fbSuppressMatching = false; // Never the less, do matching, because there are no matches
+ fbSuppressReMatching = false; // Ensure, that matching is done, if matches are absent
fTrdClusterMatches = new TClonesArray("CbmMatch", 100);
ioman->Register("TrdClusterMatch", "TRD", fTrdClusterMatches, IsOutputBranchPersistent("TrdClusterMatch"));
}
@@ -290,7 +290,7 @@ void CbmMatchRecoToMC::ReadAndCreateDataBranches()
if (nullptr != fTrdHits) {
fTrdHitMatches = static_cast<TClonesArray*>(ioman->GetObject("TrdHitMatch"));
if (nullptr == fTrdHitMatches) {
- fbSuppressMatching = false; // Never the less, do matching, because there are no matches
+ fbSuppressReMatching = false; // Ensure, that matching is done, if matches are absent
fTrdHitMatches = new TClonesArray("CbmMatch", 100);
ioman->Register("TrdHitMatch", "TRD", fTrdHitMatches, IsOutputBranchPersistent("TrdHitMatch"));
}
@@ -314,7 +314,7 @@ void CbmMatchRecoToMC::ReadAndCreateDataBranches()
if (nullptr != fMuchClusters) {
fMuchClusterMatches = static_cast<TClonesArray*>(ioman->GetObject("MuchClusterMatch"));
if (nullptr == fMuchClusterMatches) {
- fbSuppressMatching = false; // Never the less, do matching, because there are no matches
+ fbSuppressReMatching = false; // Ensure, that matching is done, if matches are absent
fMuchClusterMatches = new TClonesArray("CbmMatch", 100);
ioman->Register("MuchClusterMatch", "MUCH", fMuchClusterMatches, IsOutputBranchPersistent("MuchClusterMatch"));
}
@@ -324,7 +324,7 @@ void CbmMatchRecoToMC::ReadAndCreateDataBranches()
if (nullptr != fMuchPixelHits) {
fMuchPixelHitMatches = static_cast<TClonesArray*>(ioman->GetObject("MuchPixelHitMatch"));
if (nullptr == fMuchPixelHitMatches) {
- fbSuppressMatching = false; // Never the less, do matching, because there are no matches
+ fbSuppressReMatching = false; // Ensure, that matching is done, if matches are absent
fMuchPixelHitMatches = new TClonesArray("CbmMatch", 100);
ioman->Register("MuchPixelHitMatch", "MUCH", fMuchPixelHitMatches, IsOutputBranchPersistent("MuchPixelHitMatch"));
}
@@ -346,6 +346,7 @@ void CbmMatchRecoToMC::ReadAndCreateDataBranches()
if (nullptr != fMvdCluster) {
fMvdClusterMatches = static_cast<TClonesArray*>(ioman->GetObject("MvdClusterMatch"));
if (nullptr == fMvdClusterMatches) {
+ fbSuppressReMatching = false; // Ensure, that matching is done, if matches are absent
fMvdClusterMatches = new TClonesArray("CbmMatch", 100);
ioman->Register("MvdClusterMatch", "MVD", fMvdClusterMatches, IsOutputBranchPersistent("MvdClusterMatch"));
}
@@ -355,6 +356,7 @@ void CbmMatchRecoToMC::ReadAndCreateDataBranches()
if (nullptr != fMvdHits) {
fMvdHitMatches = static_cast<TClonesArray*>(ioman->GetObject("MvdHitMatch"));
if (nullptr == fMvdHitMatches) {
+ fbSuppressReMatching = false; // Ensure, that matching is done, if matches are absent
fMvdHitMatches = new TClonesArray("CbmMatch", 100);
ioman->Register("MvdHitMatch", "MVD", fMvdHitMatches, IsOutputBranchPersistent("MvdHitMatch"));
}
diff --git a/core/base/CbmMatchRecoToMC.h b/core/base/CbmMatchRecoToMC.h
index f093b98ec2..b686757696 100644
--- a/core/base/CbmMatchRecoToMC.h
+++ b/core/base/CbmMatchRecoToMC.h
@@ -135,14 +135,14 @@ public:
** of the branch is absent (and the corresponding cluster and hit branches are present), the flag will be set to
** false, and the cluster/hit matching will be executed.
**/
- void SuppressReMatching() { fbSuppressMatching = true; }
+ void SuppressReMatching() { fbSuppressReMatching = true; }
private:
static Int_t fEventNumber;
- Bool_t fIsMvdActive = kTRUE; // is the Mvd module active
- Bool_t fbDigiExpUsed = kTRUE; // Usage of CbmTofDigiExp instead of CbmTofDigi
- bool fbSuppressMatching = false; // Suppression of MC->hit matching, if the matches are already there
+ Bool_t fIsMvdActive = kTRUE; // is the Mvd module active
+ Bool_t fbDigiExpUsed = kTRUE; // Usage of CbmTofDigiExp instead of CbmTofDigi
+ bool fbSuppressReMatching = false; // Suppression of MC->hit matching, if the matches are already there
CbmMCDataArray* fMCTracks = nullptr; //! Monte-Carlo tracks
CbmDigiManager* fDigiManager = nullptr; //! Interface to digi branches
diff --git a/core/base/CbmTrackingDetectorInterfaceBase.h b/core/base/CbmTrackingDetectorInterfaceBase.h
index c96a2457f7..2f53e4a81c 100644
--- a/core/base/CbmTrackingDetectorInterfaceBase.h
+++ b/core/base/CbmTrackingDetectorInterfaceBase.h
@@ -148,7 +148,7 @@ public:
/// Checks detector interface: boundary conditions of the parameters
bool Check() const;
- /// Prints all the parameters into table and saves the table as a strinG
+ /// Prints all the parameters into table and saves the table as a string
std::string ToString() const;
protected:
diff --git a/core/data/tof/CbmTofHit.cxx b/core/data/tof/CbmTofHit.cxx
index 1d8cf94803..fda9f18278 100644
--- a/core/data/tof/CbmTofHit.cxx
+++ b/core/data/tof/CbmTofHit.cxx
@@ -69,6 +69,7 @@ string CbmTofHit::ToString() const
ss << "CbmTofHit: address=" << GetAddress() << " pos=(" << GetX() << "," << GetY() << "," << GetZ() << ") err=("
<< GetDx() << "," << GetDy() << "," << GetDz() << ") dxy=" << GetDxy() //<< " refId=" << GetRefId()
<< Form(" time=%8.2f", GetTime()) << " flag=" << GetFlag();
+ // << " channel=" << GetCh(); // << endl;
return ss.str();
}
diff --git a/macro/CMakeLists.txt b/macro/CMakeLists.txt
index 4cb8c24a73..32422c037a 100644
--- a/macro/CMakeLists.txt
+++ b/macro/CMakeLists.txt
@@ -2,6 +2,7 @@
add_subdirectory(run)
add_subdirectory(qa)
add_subdirectory(mcbm)
+add_subdirectory(L1)
add_subdirectory(mvd)
add_subdirectory(much)
add_subdirectory(include)
diff --git a/macro/L1/CMakeLists.txt b/macro/L1/CMakeLists.txt
new file mode 100644
index 0000000000..bbc2a4ac08
--- /dev/null
+++ b/macro/L1/CMakeLists.txt
@@ -0,0 +1 @@
+add_subdirectory(configs)
diff --git a/macro/L1/configs/CMakeLists.txt b/macro/L1/configs/CMakeLists.txt
new file mode 100644
index 0000000000..2d297127cd
--- /dev/null
+++ b/macro/L1/configs/CMakeLists.txt
@@ -0,0 +1,4 @@
+Install(FILES config_ideal_hits_mcbm.yaml
+ DESTINATION share/cbmroot/macro/L1/configs
+ )
+#Install(DIRECTORY modules DESTINATION share/cbmroot/macro/L1/configs)
diff --git a/macro/L1/configs/README b/macro/L1/configs/README
new file mode 100644
index 0000000000..94b5e815e9
--- /dev/null
+++ b/macro/L1/configs/README
@@ -0,0 +1,8 @@
+DIRECTORY: ${VMCWORKDIR}/macro/L1/configs
+
+DESCRIPTION:
+ Directory to store example configuration files, related to CA tracking.
+
+FILES:
+ config_ideal_hits_mcbm.yaml
+ Example of a configuration file for the cbm::ca::IdealHitProducer task for mCBM.
diff --git a/macro/L1/configs/config_ideal_hits_mcbm.yaml b/macro/L1/configs/config_ideal_hits_mcbm.yaml
new file mode 100644
index 0000000000..6d21e28eff
--- /dev/null
+++ b/macro/L1/configs/config_ideal_hits_mcbm.yaml
@@ -0,0 +1,42 @@
+# Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+# SPDX-License-Identifier: GPL-3.0-only
+# Authors: Sergei Zharko [committer]
+#
+## @file CbmCaIdealHitProducerConfig.yaml
+## @brief Configuration file for ideal hit producers
+## @since 28.06.2023
+## @author Sergei Zharko <s.zharko@gsi.de>
+
+ideal_hit_producer:
+ #
+ # Flags are defined for each station
+ # 0 - real hits
+ # 1 - real hit position and time are adjusted with MC (exact quantity values are used)
+ # 2 - real hit position and time are adjusted with MC (true values are smeared within corresponding error)
+ # 3 - ideal hits are created (exact quantity values are used)
+ # 4 - ideal hits are created (true values are smeared withing errors, defined in section parameters)
+ flags:
+ sts:
+ 0: 0
+ 1: 0
+ trd:
+ 0: 0
+ 1: 0
+ 2: 0
+ #
+ # Notes:
+ # 1) For MVD, MuCh, TRD and TOF du corresponds to dx, and dv corresponds to dy
+ # 2) Spatial errors are expressed in cm, time errors are expressed in ns
+ # 3) pdf stands for probability density function
+ #
+ # Available PDFs:
+ # g: bounded gaussian distribution
+ # u: uniform distribution
+ parameters:
+ sts:
+ 0: { du: { value: 1.e-3, pdf: g }, dv: { value: 1.e-3, pdf: g }, dt: { value: 5., pdf: g }}
+ 1: { du: { value: 1.e-3, pdf: g }, dv: { value: 1.e-3, pdf: g }, dt: { value: 5., pdf: g }}
+ trd:
+ 0: { du: { value: .1, pdf: g }, dv: { value: .1, pdf: g }, dt: { value: 10., pdf: g }}
+ 1: { du: { value: .1, pdf: g }, dv: { value: .1, pdf: u }, dt: { value: 10., pdf: g }}
+ 2: { du: { value: .1, pdf: u }, dv: { value: .1, pdf: g }, dt: { value: 10., pdf: g }}
diff --git a/macro/mcbm/mcbm_reco_event.C b/macro/mcbm/mcbm_reco_event.C
index b2ad563f1f..6d80f6d61a 100644
--- a/macro/mcbm/mcbm_reco_event.C
+++ b/macro/mcbm/mcbm_reco_event.C
@@ -372,22 +372,7 @@ void mcbm_reco_event(Int_t nEvents = 10, TString dataset = "data/test",
if (debugWithMC) {
auto* pIdealHitProducer = new cbm::ca::IdealHitProducer("IdealHitProducer", 3);
- //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kSts, 0, true);
- //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kSts, 1, true);
- //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kTrd, 0, false);
- //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kTrd, 1, false);
- //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kTrd, 2, false);
- //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kTof, 0, false);
- //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kTof, 1, false);
- //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kTof, 2, false);
- //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kTof, 3, false);
- //pIdealHitProducer->CreateNewHits(L1DetectorID::kSts, 0, true, 1.e-3, 1.e-3, 5.);
- //pIdealHitProducer->CreateNewHits(L1DetectorID::kSts, 1, true, 1.e-3, 1.e-3, 5.);
- //pIdealHitProducer->CreateNewHits(L1DetectorID::kTrd, 0, true, .1, .1, 10.);
- //pIdealHitProducer->CreateNewHits(L1DetectorID::kTrd, 1, true, .1, .1, 10.);
- //pIdealHitProducer->CreateNewHits(L1DetectorID::kTrd, 2, false, .1, .1, 10.);
- //pIdealHitProducer->CreateNewHits(L1DetectorID::kTrd, 1, false, .1, .1, 10.);
- //pIdealHitProducer->CreateNewHits(L1DetectorID::kTrd, 2, false, .1, .1, 10.);
+ pIdealHitProducer->SetConfigName(srcDir + "/macro/L1/configs/config_ideal_hits_mcbm.yaml");
//run->AddTask(pIdealHitProducer);
}
diff --git a/macro/qa/qa_compare_ca.C b/macro/qa/qa_compare_ca.C
new file mode 100644
index 0000000000..3bc59effd5
--- /dev/null
+++ b/macro/qa/qa_compare_ca.C
@@ -0,0 +1,45 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// @file qa_compare_ca.C
+/// @author Sergei Zharko <s.zharko@gsi.de>
+/// @since 22.06.2023
+/// @brief ROOT macro to run QA-Checker framework comparison of the old and new ROOT-file versions (CA).
+
+/// @brief Function to compare two different versions of ROOT-files
+/// @param configName Name of config containing file-object map
+/// @param outputName Name of the QA-Checker output (default is "QaCheckerResult.root")
+/// @return Result flag:
+/// - 0: all objects are the same within defined comparison procedure (point-by-point for now)
+/// - 1: some of the objects differ
+
+
+/// NOTE: Disable clang formatting to keep easy parameters overview
+/* clang-format off */
+int qa_compare_ca(
+ const char* configName = "objects2.yaml",
+ const char* outputName = "QACheckerOutput.root" // TODO: Add tag of a merge request or commit
+ )
+/* clang-format on */
+{
+ // ----- Logger settings
+ FairLogger::GetLogger()->SetLogScreenLevel("INFO");
+ FairLogger::GetLogger()->SetColoredLog(true);
+
+ //// ----- Style settings
+ gStyle->SetPalette(kRainBow);
+
+ //// ----- Configure QA-Checker
+ auto pQaChecker = std::make_unique<cbm::qa::checker::Core>();
+ pQaChecker->RegisterOutFile(outputName); // Set name of the output file
+ pQaChecker->SetFromYAML(configName); // Read file-object map
+
+ //// ----- Run comparision routine
+ pQaChecker->Process("P");
+
+ //// ----- Scan results
+ bool res = pQaChecker->Scan(); // true - objects are the same, false - objects differ
+ std::cout << "Macro finished successfully." << std::endl;
+ return !res;
+}
diff --git a/macro/run/run_reco.C b/macro/run/run_reco.C
index 71b5c867fb..9f3be86520 100644
--- a/macro/run/run_reco.C
+++ b/macro/run/run_reco.C
@@ -387,40 +387,6 @@ void run_reco(TString input = "", Int_t nTimeSlices = -1, Int_t firstTimeSlice =
auto kalman = new CbmKF();
run->AddTask(kalman);
- if (debugWithMC) {
- auto* pIdealHitProducer = new cbm::ca::IdealHitProducer("IdealHitProducer", 3);
- for (int iSt = 0; iSt < 8; ++iSt) {
- if (iSt == 2) continue;
- pIdealHitProducer->CreateNewHits(L1DetectorID::kSts, iSt, false, 1.e-3, 1.e-3, 5.);
- }
- //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kSts, 0, false);
- //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kSts, 1, false);
- //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kSts, 2, false);
- //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kSts, 3, false);
- //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kSts, 4, false);
- //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kSts, 5, false);
- //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kSts, 6, false);
- //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kSts, 7, false);
- //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kMuch, 0, false);
- //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kMuch, 1, false);
- //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kMuch, 2, false);
- //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kMuch, 3, false);
- //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kMuch, 4, false);
- //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kMuch, 5, false);
- //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kMuch, 6, false);
- //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kMuch, 7, false);
- //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kMuch, 8, false);
- //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kMuch, 9, false);
- //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kMuch, 10, false);
- //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kMuch, 11, false);
- //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kTrd, 0, false);
- //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kTrd, 1, false);
- //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kTrd, 2, false);
- //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kTrd, 3, false);
- //pIdealHitProducer->AdjustExistingHits(L1DetectorID::kTof, 0, false);
- //run->AddTask(pIdealHitProducer);
- }
-
// L1 tracking
auto l1 = (debugWithMC) ? new CbmL1("CA", 2, 3) : new CbmL1("CA");
diff --git a/reco/L1/L1Algo/utils/CaAlgoRandom.h b/reco/L1/L1Algo/utils/CaAlgoRandom.h
index d506496c5c..ba3dd2a702 100644
--- a/reco/L1/L1Algo/utils/CaAlgoRandom.h
+++ b/reco/L1/L1Algo/utils/CaAlgoRandom.h
@@ -59,13 +59,21 @@ namespace ca::algo
void SetSeed(int seed);
/// @brief Returns a normally distributed random value, limited within a selected sigma range
- /// @tparam T Type of floating point
+ /// @tparam T Type of floating point numbers
/// @param mean Mean of the distribution
/// @param sigma Sigma of the distribution
/// @param nSigmas Half-width of the generated numbers domain, expressed in number of sigmas
template<typename T, std::enable_if_t<std::is_floating_point<T>::value, bool> = true>
T BoundedGaus(const T& mean, const T& sigma, const T& nSigmas) const;
+ /// @brief Returns a random value, addressed with a continuous uniform distribution
+ /// @tparam T Type of floating point numbers
+ /// @param mean Mean of the distribution
+ /// @param sigma Sigma of the distribution
+ template<typename T, std::enable_if_t<std::is_floating_point<T>::value, bool> = true>
+ T Uniform(const T& mean, const T& sigma) const;
+
+
private:
mutable GeneratorType_t fGenerator; ///< Random number generator
unsigned int fSeed = 1; ///< Random number seed
@@ -93,5 +101,15 @@ T ca::algo::Random::BoundedGaus(const T& mean, const T& sigma, const T& nSigmas)
return res;
}
+// ---------------------------------------------------------------------------------------------------------------------
+//
+template<typename T, std::enable_if_t<std::is_floating_point<T>::value, bool>>
+T ca::algo::Random::Uniform(const T& mean, const T& sigma) const
+{
+ LOG_IF(fatal, !(sigma > 0 && std::isfinite(sigma))) << "ca::algo::Random::Uniform sigma = " << sigma << " is illegal";
+ std::uniform_real_distribution rnd {-sigma * std::sqrt(3) + mean, sigma * std::sqrt(3) + mean};
+ return rnd(fGenerator);
+}
+
#endif // CaAlgoRandom_h
diff --git a/reco/L1/utils/CbmCaIdealHitProducer.cxx b/reco/L1/utils/CbmCaIdealHitProducer.cxx
index 873787954d..441dc328b2 100644
--- a/reco/L1/utils/CbmCaIdealHitProducer.cxx
+++ b/reco/L1/utils/CbmCaIdealHitProducer.cxx
@@ -11,40 +11,11 @@
#include "CbmSetup.h"
+
using cbm::ca::IdealHitProducer;
ClassImp(cbm::ca::IdealHitProducer);
-// ---------------------------------------------------------------------------------------------------------------------
-//
-void IdealHitProducer::AdjustExistingHits(L1DetectorID detID, int iSt, bool ifSmear)
-{
- switch (detID) {
- case L1DetectorID::kMvd: fHitProducerMvd.AdjustExistingHits(iSt, ifSmear); break;
- case L1DetectorID::kSts: fHitProducerSts.AdjustExistingHits(iSt, ifSmear); break;
- case L1DetectorID::kMuch: fHitProducerMuch.AdjustExistingHits(iSt, ifSmear); break;
- case L1DetectorID::kTrd: fHitProducerTrd.AdjustExistingHits(iSt, ifSmear); break;
- case L1DetectorID::kTof: fHitProducerTof.AdjustExistingHits(iSt, ifSmear); break;
- case L1DetectorID::kEND:
- LOG(fatal) << fName << "::AdjustExistingHits(): L1DetectorID::kEND cannot be passed as a detID";
- }
-}
-
-// ---------------------------------------------------------------------------------------------------------------------
-//
-void IdealHitProducer::CreateNewHits(L1DetectorID detID, int iSt, bool ifSmear, double du, double dv, double dt)
-{
- switch (detID) {
- case L1DetectorID::kMvd: fHitProducerMvd.CreateNewHits(iSt, ifSmear, du, dv, dt); break;
- case L1DetectorID::kSts: fHitProducerSts.CreateNewHits(iSt, ifSmear, du, dv, dt); break;
- case L1DetectorID::kMuch: fHitProducerMuch.CreateNewHits(iSt, ifSmear, du, dv, dt); break;
- case L1DetectorID::kTrd: fHitProducerTrd.CreateNewHits(iSt, ifSmear, du, dv, dt); break;
- case L1DetectorID::kTof: fHitProducerTof.CreateNewHits(iSt, ifSmear, du, dv, dt); break;
- case L1DetectorID::kEND:
- LOG(fatal) << fName << "::AdjustExistingHits(): L1DetectorID::kEND cannot be passed as a detID";
- }
-}
-
// ---------------------------------------------------------------------------------------------------------------------
//
InitStatus IdealHitProducer::Init()
@@ -75,3 +46,14 @@ void IdealHitProducer::Exec(Option_t* option)
if (fbUseDet[L1DetectorID::kTrd]) { fHitProducerTrd.Exec(option); }
if (fbUseDet[L1DetectorID::kTof]) { fHitProducerTof.Exec(option); }
}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void IdealHitProducer::SetConfigName(const char* name)
+{
+ fHitProducerMvd.SetConfigName(name);
+ fHitProducerSts.SetConfigName(name);
+ fHitProducerMuch.SetConfigName(name);
+ fHitProducerTrd.SetConfigName(name);
+ fHitProducerTof.SetConfigName(name);
+}
diff --git a/reco/L1/utils/CbmCaIdealHitProducer.h b/reco/L1/utils/CbmCaIdealHitProducer.h
index 6530c0de5a..0b5e175317 100644
--- a/reco/L1/utils/CbmCaIdealHitProducer.h
+++ b/reco/L1/utils/CbmCaIdealHitProducer.h
@@ -15,6 +15,8 @@
#include "FairTask.h"
+#include <string>
+
namespace cbm::ca
{
/// @brief Ideal hit producer task for CA tracking
@@ -40,21 +42,6 @@ namespace cbm::ca
/// @brief Move assignment operator
IdealHitProducer& operator=(IdealHitProducer&&) = delete;
- /// @brief Sets a mode to adjust existing hits to MC points for a selected tracking station
- /// @param detID Detector ID
- /// @param iSt Index of station
- /// @param ifSmear Flag: true - MC true position and time are smeared within error
- void AdjustExistingHits(L1DetectorID detID, int iSt, bool ifSmear);
-
- /// @brief Sets a mode to create hits from MC points for a selected tracking station
- /// @param detID Detector ID
- /// @param iSt Index of station
- /// @param ifSmear Flag: true - MC true position and time are smeared within error
- /// @param du Error of u-coordinate measurement [cm]
- /// @param dv Error of v-coordinate measurement [cm]
- /// @param dt Error of time measurement [ns]
- void CreateNewHits(L1DetectorID detID, int iSt, bool ifSmear, double du, double dv, double dt);
-
/// @brief Initialization of the task
InitStatus Init();
@@ -64,6 +51,10 @@ namespace cbm::ca
/// @brief Execution of the task
void Exec(Option_t* option);
+ /// @brief Sets YAML configuration file with defined smearing parameters
+ /// @param name Name of the configuration file
+ void SetConfigName(const char* name);
+
ClassDef(IdealHitProducer, 1);
private:
diff --git a/reco/L1/utils/CbmCaIdealHitProducerDet.h b/reco/L1/utils/CbmCaIdealHitProducerDet.h
index cd40487f4f..b4f28a8be3 100644
--- a/reco/L1/utils/CbmCaIdealHitProducerDet.h
+++ b/reco/L1/utils/CbmCaIdealHitProducerDet.h
@@ -36,10 +36,7 @@
#include "CbmTrdPoint.h"
#include "CbmTrdTrackingInterface.h"
-<<<<<<< HEAD
#include "FairRootManager.h"
-=======
->>>>>>> CA: Bugfix in ideal hit finder
#include "FairTask.h"
#include "Logger.h"
@@ -47,9 +44,14 @@
#include "TVector3.h"
#include <array>
+#include <optional>
+#include <sstream>
+#include <string>
#include <tuple>
#include <unordered_map>
+#include <yaml-cpp/yaml.h>
+
#include "CaAlgoRandom.h"
#include "L1Constants.h"
#include "L1Def.h"
@@ -84,19 +86,6 @@ namespace cbm::ca
/// @brief Move assignment operator
IdealHitProducerDet& operator=(IdealHitProducerDet&&) = delete;
- /// @brief Sets a mode to adjust existing hits to MC points for a selected tracking station
- /// @param iSt Index of station
- /// @param ifSmear Flag: true - MC true position and time are smeared within error
- void AdjustExistingHits(int iSt, bool ifSmear);
-
- /// @brief Sets a mode to create hits from MC points for a selected tracking station
- /// @param iSt Index of station
- /// @param ifSmear Flag: true - MC true position and time are smeared within error
- /// @param du Error of u-coordinate measurement [cm]
- /// @param dv Error of v-coordinate measurement [cm]
- /// @param dt Error of time measurement [ns]
- void CreateNewHits(int iSt, bool ifSmear, double du, double dv, double dt);
-
/// @brief Initialization of the task
InitStatus Init();
@@ -110,6 +99,10 @@ namespace cbm::ca
/// @param seed Random seed
void SetRandomSeed(int seed) { fRandom.SetSeed(seed); }
+ /// @brief Sets YAML configuration file name
+ /// @param name Name of the configuration file
+ void SetConfigName(const char* name) { fsConfigName = name; }
+
private:
static constexpr double kHitWSize = 3.5; ///< half-width for the bounded gaussian [sigma]
@@ -118,34 +111,42 @@ namespace cbm::ca
/// @return link to matched point
std::optional<CbmLink> GetMatchedPointLink(int iH);
+ /// @brief Parses the YAML configuration file
+ void ParseConfig();
+
/// @brief Pushes back a hit into the hits branch
/// @param pHit Pointer to source hit
void PushBackHit(const Hit_t* pHit);
+ /// @brief Smears the value by a given standard deviation with a selected function
+ /// @param value Value to be smeared
+ /// @param sigma Value of the standard deviation
+ /// @param opt Function (0 - BoundedGaus, 1 - Uniform)
+ void SmearValue(double& value, double sigma, int opt);
/// @brief A structure to keep hit adjustment/creation settings for each station
- struct StationSetting {
+ struct HitParameters {
double fDu = undef::kD; ///< Error of u-coordinate measurement [cm]
double fDv = undef::kD; ///< Error of v-coordinate measurement [cm]
double fDt = undef::kD; ///< Error of time measurement [ns]
+ /// PDFs selection for different quantities
+ /// - 0: Bounded Gaussian distribution
+ /// - 1: Uniform distribution
+ short fPdfU = -1;
+ short fPdfV = -1;
+ short fPdfT = -1;
+
/// @brief Mode:
- /// - 0: do nothing with hits
- /// - 1: adjust measured hits with matched point info (hit errors used)
- /// - 2: create new hits from points
+ /// - 0: real hits
+ /// - 1: existing hits are improved from MC points
+ /// - 2: new hits are created from MC points
short fMode = 0;
/// @brief Flag:
/// - true: MC point quantities are smeared with gaussian using error
/// - false: Precise MC point quantities are used
bool fbSmear = false;
-
- /// @brief String representation of the structure object
- /// @param verbose Verbosity level:
- /// - 0: prints nothing
- /// - 1: prints one-line log
- /// - 2: prints verbose log
- std::string ToString(int verbose = 1) const;
};
CbmTrackingDetectorInterfaceBase* fpDetInterface = nullptr; ///< Detector interface
@@ -162,13 +163,16 @@ namespace cbm::ca
TClonesArray* fpBrHitsTmp = nullptr; ///< Temporary array of hits
TClonesArray* fpBrHitMatchesTmp = nullptr; ///< Temporary array of hit matches
- std::unordered_map<int, StationSetting> fmStationPars; ///< List of settings vs. tracking station
+ std::string fsConfigName = ""; ///< Name of configuration file
+
+ std::vector<HitParameters> fvStationPars; ///< Parameters, stored for each station
::ca::algo::Random fRandom {1}; ///< Random generator
/// @brief Management flag, which does run the routine if there was a request
bool fbRunTheRoutine = true;
+
int fHitCounter = 0; ///< Hit counter in a new branch
};
@@ -195,31 +199,6 @@ namespace cbm::ca
}
}
- // ------------------------------------------------------------------------------------------------------------------
- //
- template<L1DetectorID DetID>
- void IdealHitProducerDet<DetID>::AdjustExistingHits(int iSt, bool ifSmear)
- {
- StationSetting pars {};
- pars.fMode = 1;
- pars.fbSmear = ifSmear;
- fmStationPars[iSt] = pars;
- }
-
- // ------------------------------------------------------------------------------------------------------------------
- //
- template<L1DetectorID DetID>
- void IdealHitProducerDet<DetID>::CreateNewHits(int iSt, bool ifSmear, double du, double dv, double dt)
- {
- StationSetting pars {};
- pars.fMode = 2;
- pars.fbSmear = ifSmear;
- pars.fDu = du;
- pars.fDv = dv;
- pars.fDt = dt;
- fmStationPars[iSt] = pars;
- }
-
// ------------------------------------------------------------------------------------------------------------------
//
template<L1DetectorID DetID>
@@ -293,26 +272,21 @@ namespace cbm::ca
}
// ------ Checks of the settings of hit creation/adjustment mode
- int nStations = fpDetInterface->GetNtrackingStations();
- for (int iSt = 0; iSt < nStations; ++iSt) {
- if (fmStationPars.find(iSt) == fmStationPars.end()) {
- StationSetting pars {};
- pars.fMode = 0;
- fmStationPars[iSt] = pars;
- }
+ if (fsConfigName.size() != 0) {
+ this->ParseConfig();
+ // Skip event/time slice processing, if there are no requests for hit modification of the detector
+ // (fbRunTheRoutine = true, if there is at least one station with fMode > 0)
+ fbRunTheRoutine =
+ std::any_of(fvStationPars.begin(), fvStationPars.end(), [](const auto& p) { return p.fMode > 0; });
+ }
+ else {
+ fbRunTheRoutine = false;
}
-
- // Skip event/time slice processing, if there are no requests for hit modification of the detector
- // (fbRunTheRoutine = true, if there is at least one station with fMode > 0)
- fbRunTheRoutine =
- std::any_of(fmStationPars.begin(), fmStationPars.end(), [](const auto& p) { return p.second.fMode > 0; });
std::stringstream msg;
msg << "\033[1;31mATTENTION! IDEAL HIT PRODUCER IS USED FOR " << kDetName[DetID] << "\033[0m";
-
LOG_IF(info, fbRunTheRoutine) << msg.str();
-
return kSUCCESS;
}
catch (const std::logic_error& error) {
@@ -321,6 +295,94 @@ namespace cbm::ca
return kERROR;
}
+ // -------------------------------------------------------------------------------------------------------------------
+ //
+ template<L1DetectorID DetID>
+ void IdealHitProducerDet<DetID>::ParseConfig()
+ {
+ int nStations = fpDetInterface->GetNtrackingStations();
+ fvStationPars.clear();
+ fvStationPars.resize(nStations);
+
+ // Read file
+ YAML::Node config;
+ try {
+ config = YAML::LoadFile(fsConfigName)["ideal_hit_producer"];
+
+ std::string detBranchName = "";
+ if constexpr (L1DetectorID::kMvd == DetID) { detBranchName = "mvd"; }
+ else if constexpr (L1DetectorID::kSts == DetID) {
+ detBranchName = "sts";
+ }
+ else if constexpr (L1DetectorID::kMuch == DetID) {
+ detBranchName = "much";
+ }
+ else if constexpr (L1DetectorID::kTrd == DetID) {
+ detBranchName = "trd";
+ }
+ else if constexpr (L1DetectorID::kTof == DetID) {
+ detBranchName = "tof";
+ }
+
+ const auto& parNode = config["parameters"][detBranchName];
+ const auto& flagNode = config["flags"][detBranchName];
+
+ for (int iSt = 0; iSt < nStations; ++iSt) {
+ auto& par = fvStationPars[iSt];
+ auto& node = parNode[iSt];
+ short mode = flagNode[iSt].as<int>(0);
+ if (mode > 0) {
+ par.fDu = node["du"]["value"].as<double>();
+ switch (node["du"]["pdf"].as<char>()) {
+ case 'g': par.fPdfU = 0; break;
+ case 'u': par.fPdfU = 1; break;
+ default: par.fPdfU = -1; break;
+ }
+ par.fDv = node["dv"]["value"].as<double>();
+ switch (node["dv"]["pdf"].as<char>()) {
+ case 'g': par.fPdfV = 0; break;
+ case 'u': par.fPdfV = 1; break;
+ default: par.fPdfV = -1; break;
+ }
+ par.fDt = node["dt"]["value"].as<double>();
+ switch (node["dt"]["pdf"].as<char>()) {
+ case 'g': par.fPdfT = 0; break;
+ case 'u': par.fPdfT = 1; break;
+ default: par.fPdfT = -1; break;
+ }
+ switch (mode) {
+ case 1:
+ par.fMode = 1;
+ par.fbSmear = false;
+ break;
+ case 2:
+ par.fMode = 1;
+ par.fbSmear = true;
+ break;
+ case 3:
+ par.fMode = 2;
+ par.fbSmear = false;
+ break;
+ case 4:
+ par.fMode = 2;
+ par.fbSmear = true;
+ break;
+ }
+ }
+ }
+ }
+ catch (const YAML::BadFile& err) {
+ std::stringstream msg;
+ msg << "YAML configuration file " << fsConfigName << " does not exist";
+ throw std::runtime_error(msg.str());
+ }
+ catch (const YAML::ParserException& err) {
+ std::stringstream msg;
+ msg << "YAML configuration file " << fsConfigName << " is formatted improperly";
+ throw std::runtime_error(msg.str());
+ }
+ }
+
// -------------------------------------------------------------------------------------------------------------------
//
template<L1DetectorID DetID>
@@ -393,6 +455,18 @@ namespace cbm::ca
return std::nullopt;
}
+ // -------------------------------------------------------------------------------------------------------------------
+ //
+ template<L1DetectorID DetID>
+ void IdealHitProducerDet<DetID>::SmearValue(double& value, double sigma, int opt)
+ {
+ switch (opt) {
+ case 0: value += fRandom.BoundedGaus<double>(0., sigma, kHitWSize); break;
+ case 1: value += fRandom.Uniform<double>(0., sigma); break;
+ default: LOG(fatal) << "IdealHitProducerDet::SmearValue: illegal option: opt = " << opt;
+ }
+ }
+
// -------------------------------------------------------------------------------------------------------------------
//
template<L1DetectorID DetID>
@@ -426,7 +500,7 @@ namespace cbm::ca
// Get setting
int iSt = fpDetInterface->GetTrackingStationIndex(pHit);
- const auto& setting = fmStationPars.at(iSt);
+ const auto& setting = fvStationPars[iSt];
// Skip hits from stations, where new hits are to be created from points
if (setting.fMode == 2) { continue; }
@@ -462,17 +536,18 @@ namespace cbm::ca
double t = pPoint->GetTime() + fpMCEventList->GetEventTime(link);
if (setting.fbSmear) {
- auto dt = pHit->GetTimeError();
- auto du = pHit->GetDx();
- auto dv = pHit->GetDy();
+ double dt = pHit->GetTimeError();
+ double du = pHit->GetDx();
+ double dv = pHit->GetDy();
if constexpr (L1DetectorID::kSts == DetID) {
du = pHit->GetDu();
dv = pHit->GetDv();
}
auto [u, v] = fpDetInterface->ConvertXYtoUV(iSt, x, y);
- u += fRandom.BoundedGaus<double>(0., du, kHitWSize);
- v += fRandom.BoundedGaus<double>(0., dv, kHitWSize);
- t += fRandom.BoundedGaus<double>(0., dt, kHitWSize);
+ this->SmearValue(u, du, setting.fPdfU);
+ this->SmearValue(v, dv, setting.fPdfV);
+ this->SmearValue(t, dt, setting.fPdfT);
+
std::tie(x, y) = fpDetInterface->ConvertUVtoXY(iSt, u, v);
}
pHit->SetX(x);
@@ -505,7 +580,7 @@ namespace cbm::ca
for (int iP = 0; iP < nPoints; ++iP) {
auto* pPoint = static_cast<Point_t*>(fpBrPoints->Get(fileId, eventId, iP));
int iSt = fpDetInterface->GetTrackingStationIndex(pPoint->GetDetectorID());
- const auto& setting = fmStationPars.at(iSt);
+ const auto& setting = fvStationPars[iSt];
if (setting.fMode != 2) { continue; }
double du = setting.fDu;
double dv = setting.fDv;
@@ -534,9 +609,9 @@ namespace cbm::ca
if (setting.fbSmear) {
auto [u, v] = fpDetInterface->ConvertXYtoUV(iSt, x, y);
// TODO: Provide more realistic profiles for TRD
- u += fRandom.BoundedGaus<double>(0., du, kHitWSize);
- v += fRandom.BoundedGaus<double>(0., dv, kHitWSize);
- t += fRandom.BoundedGaus<double>(0., dt, kHitWSize);
+ this->SmearValue(u, du, setting.fPdfU);
+ this->SmearValue(v, dv, setting.fPdfV);
+ this->SmearValue(t, dt, setting.fPdfT);
std::tie(x, y) = fpDetInterface->ConvertUVtoXY(iSt, u, v);
}
pos.SetX(x);
--
GitLab