diff --git a/core/qa/CbmQaCmpDrawer.cxx b/core/qa/CbmQaCmpDrawer.cxx
index 84ef5bc25b80c547b79e812f32e6fe6c31548b5a..ed2b2305db77f7a965137be7ee3b3d1765862d69 100644
--- a/core/qa/CbmQaCmpDrawer.cxx
+++ b/core/qa/CbmQaCmpDrawer.cxx
@@ -18,4 +18,3 @@ template class CbmQaCmpDrawer<TH1F>;
template class CbmQaCmpDrawer<TProfile>;
template class CbmQaCmpDrawer<TH2F>;
template class CbmQaCmpDrawer<TProfile2D>;
-
diff --git a/macro/mcbm/mcbm_qa.C b/macro/mcbm/mcbm_qa.C
index 363d9c0059d99f2c3d18b8cdb34cceedb5111523..2d4d45fd99a7148c8910f03cded9b10348475f6d 100644
--- a/macro/mcbm/mcbm_qa.C
+++ b/macro/mcbm/mcbm_qa.C
@@ -244,7 +244,7 @@ void mcbm_qa(Int_t nEvents = 0, TString dataset = "data/mcbm_beam_2020_03_test",
pCaOutputQa->SetMcbmTrackingMode();
pCaOutputQa->ReadParameters(caParFile.Data());
pCaOutputQa->SetUseSts(bUseSts);
- //pCaOutputQa->SetUseMuch(bUseMuch);
+ pCaOutputQa->SetUseMuch(bUseMuch);
pCaOutputQa->SetUseTrd(bUseTrd);
pCaOutputQa->SetUseTof(bUseTof);
run->AddTask(pCaOutputQa);
diff --git a/reco/L1/CMakeLists.txt b/reco/L1/CMakeLists.txt
index 53086f8c772ac6acef3623a216fd2b7294c177fa..469c472599a677bccae98991fa0b832735459aa0 100644
--- a/reco/L1/CMakeLists.txt
+++ b/reco/L1/CMakeLists.txt
@@ -38,7 +38,7 @@ set(SRCS
OffLineInterface/CbmL1StsTrackFinder.cxx
OffLineInterface/CbmL1GlobalTrackFinder.cxx
OffLineInterface/CbmL1GlobalFindTracksEvents.cxx
-
+ L1Algo/utils/CaAlgoRandom.cxx
L1Algo/L1Algo.cxx
L1Algo/L1TrackPar.cxx
L1Algo/L1CaTrackFinder.cxx
@@ -115,6 +115,7 @@ set(HEADERS
L1Algo/L1Undef.h
catools/CaToolsWindowFinder.h
catools/CaToolsLinkKey.h
+ catools/CaToolsHitRecord.h
qa/CbmCaInputQaBase.h
)
diff --git a/reco/L1/CbmCaMCModule.cxx b/reco/L1/CbmCaMCModule.cxx
index 2263509a63b89985b74b0d7d53927c0dbe542e80..db1e8fc942f59fd4df8f048f8c97f1e3020f2b02 100644
--- a/reco/L1/CbmCaMCModule.cxx
+++ b/reco/L1/CbmCaMCModule.cxx
@@ -1,4 +1,4 @@
-/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+/* Copyright (C) 2022-2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
SPDX-License-Identifier: GPL-3.0-only
Authors: Sergei Zharko [committer] */
@@ -19,7 +19,6 @@
#include "CbmMCTrack.h"
#include "CbmStsHit.h"
#include "CbmTimeSlice.h"
-
#include "FairEventHeader.h"
#include "FairMCEventHeader.h"
#include "FairRootManager.h"
@@ -42,6 +41,8 @@
// *********************************
using cbm::ca::MCModule;
+using L1Constants::clrs::kCL; // clear log
+using L1Constants::clrs::kRDb; // red bold log
// ---------------------------------------------------------------------------------------------------------------------
//
@@ -63,60 +64,43 @@ try {
fpMCEventHeader = mcManager->GetObject("MCEventHeader.");
fpMCTracks = mcManager->InitBranch("MCTrack");
- fpMvdPoints = nullptr;
- fpStsPoints = nullptr;
- fpMuchPoints = nullptr;
- fpTrdPoints = nullptr;
- fpTofPoints = nullptr;
-
- fpMvdHitMatches = nullptr;
- fpStsHitMatches = nullptr;
- fpMuchHitMatches = nullptr;
- fpTrdHitMatches = nullptr;
- fpTofHitMatches = nullptr;
+ fvpBrPoints.fill(nullptr);
+ fvpBrHitMatches.fill(nullptr);
fpStsClusterMatches = nullptr;
fpStsHits = nullptr;
fFileEventIDs.clear();
- if (fbUseMvd) {
- LOG(info) << "CA MC Module: initializing branches for MVD";
- fpMvdPoints = mcManager->InitBranch("MvdPoint");
- fpMvdHitMatches = dynamic_cast<TClonesArray*>(fairManager->GetObject("MvdHitMatch"));
- }
+ auto InitPointBranch = [&](const char* brName, L1DetectorID detID) {
+ if (fvbUseDet[detID]) { fvpBrPoints[detID] = mcManager->InitBranch(brName); }
+ };
- if (fbUseSts) {
- LOG(info) << "CA MC Module: initializing branches for STS";
- fpStsPoints = mcManager->InitBranch("StsPoint");
- fpStsHitMatches = dynamic_cast<TClonesArray*>(fairManager->GetObject("StsHitMatch"));
- fpStsClusterMatches = dynamic_cast<TClonesArray*>(fairManager->GetObject("StsClusterMatch"));
- fpStsHits = dynamic_cast<TClonesArray*>(fairManager->GetObject("StsHit"));
- }
+ auto InitMatchesBranch = [&](const char* brName, L1DetectorID detID) {
+ if (fvbUseDet[detID]) { fvpBrHitMatches[detID] = dynamic_cast<TClonesArray*>(fairManager->GetObject(brName)); }
+ };
- if (fbUseMuch) {
- LOG(info) << "CA MC Module: initializing branches for MuCh";
- fpMuchPoints = mcManager->InitBranch("MuchPoint");
- fpMuchHitMatches = dynamic_cast<TClonesArray*>(fairManager->GetObject("MuchPixelHitMatch"));
- }
+ InitPointBranch("MvdPoint", L1DetectorID::kMvd);
+ InitPointBranch("StsPoint", L1DetectorID::kSts);
+ InitPointBranch("MuchPoint", L1DetectorID::kMuch);
+ InitPointBranch("TrdPoint", L1DetectorID::kTrd);
+ InitPointBranch("TofPoint", L1DetectorID::kTof);
- if (fbUseTrd) {
- LOG(info) << "CA MC Module: initializing branches for TRD";
- fpTrdPoints = mcManager->InitBranch("TrdPoint");
- fpTrdHitMatches = dynamic_cast<TClonesArray*>(fairManager->GetObject("TrdHitMatch"));
- }
+ InitMatchesBranch("MvdHitMatch", L1DetectorID::kMvd);
+ InitMatchesBranch("StsHitMatch", L1DetectorID::kSts);
+ InitMatchesBranch("MuchPixelHitMatch", L1DetectorID::kMuch);
+ InitMatchesBranch("TrdHitMatch", L1DetectorID::kTrd);
+ InitMatchesBranch("TofHitMatch", L1DetectorID::kTof);
- if (fbUseTof) {
- LOG(info) << "CA MC Module: initializing branches for TOF";
- fpTofPoints = mcManager->InitBranch("TofPoint");
- fpTofHitMatches = dynamic_cast<TClonesArray*>(fairManager->GetObject("TofHitMatch"));
+ if (fvbUseDet[L1DetectorID::kSts]) {
+ fpStsClusterMatches = dynamic_cast<TClonesArray*>(fairManager->GetObject("StsClusterMatch"));
+ fpStsHits = dynamic_cast<TClonesArray*>(fairManager->GetObject("StsHit"));
}
// Check initialization
this->CheckInit();
LOG(info) << "CA MC Module: initializing CA tracking Monte-Carlo module... \033[1;32mDone!\033[0m";
-
return true;
}
catch (const std::logic_error& error) {
@@ -217,7 +201,11 @@ void MCModule::Finish() {}
//
void MCModule::MatchRecoAndMC()
{
- this->MatchPointsAndHits();
+ this->MatchPointsAndHits<L1DetectorID::kMvd>();
+ this->MatchPointsAndHits<L1DetectorID::kSts>();
+ this->MatchPointsAndHits<L1DetectorID::kMuch>();
+ this->MatchPointsAndHits<L1DetectorID::kTrd>();
+ this->MatchPointsAndHits<L1DetectorID::kTof>();
this->MatchRecoAndMCTracks();
this->InitTrackInfo();
}
@@ -229,13 +217,10 @@ template<>
int MCModule::MatchHitWithMc<L1DetectorID::kMvd>(int iHitExt) const
{
int iPoint = -1;
- if (fpMvdHitMatches) {
- const auto* hitMatch = dynamic_cast<CbmMatch*>(fpMvdHitMatches->At(iHitExt));
- assert(hitMatch);
- if (hitMatch->GetNofLinks() > 0
- && hitMatch->GetLink(0).GetIndex() < fpMCData->GetNofPointsUsed(L1DetectorID::kMvd)) {
- iPoint = hitMatch->GetLink(0).GetIndex();
- }
+ const auto* hitMatch = dynamic_cast<CbmMatch*>(fvpBrHitMatches[L1DetectorID::kMvd]->At(iHitExt));
+ assert(hitMatch);
+ if (hitMatch->GetNofLinks() > 0 && hitMatch->GetLink(0).GetIndex() < fpMCData->GetNofPointsUsed(L1DetectorID::kMvd)) {
+ iPoint = hitMatch->GetLink(0).GetIndex();
}
return iPoint;
}
@@ -302,7 +287,7 @@ template<>
int MCModule::MatchHitWithMc<L1DetectorID::kMuch>(int iHitExt) const
{
int iPoint = -1;
- const auto* hitMatchMuch = dynamic_cast<CbmMatch*>(fpMuchHitMatches->At(iHitExt));
+ const auto* hitMatchMuch = dynamic_cast<CbmMatch*>(fvpBrHitMatches[L1DetectorID::kMuch]->At(iHitExt));
if (hitMatchMuch) {
for (int iLink = 0; iLink < hitMatchMuch->GetNofLinks(); ++iLink) {
if (hitMatchMuch->GetLink(iLink).GetIndex() < fpMCData->GetNofPointsUsed(L1DetectorID::kMuch)) {
@@ -327,7 +312,7 @@ template<>
int MCModule::MatchHitWithMc<L1DetectorID::kTrd>(int iHitExt) const
{
int iPoint = -1;
- const auto* hitMatch = dynamic_cast<const CbmMatch*>(fpTrdHitMatches->At(iHitExt));
+ const auto* hitMatch = dynamic_cast<const CbmMatch*>(fvpBrHitMatches[L1DetectorID::kTrd]->At(iHitExt));
if (hitMatch) {
int iMC = -1;
if (hitMatch->GetNofLinks() > 0) {
@@ -351,8 +336,9 @@ int MCModule::MatchHitWithMc<L1DetectorID::kTrd>(int iHitExt) const
template<>
int MCModule::MatchHitWithMc<L1DetectorID::kTof>(int iHitExt) const
{
+ // TODO: SZh 23.05.2023: Use TOF interactions class
int iPoint = -1;
- const auto* hitMatch = dynamic_cast<const CbmMatch*>(fpTofHitMatches->At(iHitExt));
+ const auto* hitMatch = dynamic_cast<const CbmMatch*>(fvpBrHitMatches[L1DetectorID::kTof]->At(iHitExt));
if (hitMatch) {
for (int iLink = 0; iLink < hitMatch->GetNofLinks(); ++iLink) {
int iMc = hitMatch->GetLink(iLink).GetIndex();
@@ -360,7 +346,7 @@ int MCModule::MatchHitWithMc<L1DetectorID::kTof>(int iHitExt) const
int iEvent = hitMatch->GetLink(iLink).GetEntry();
if (iMc < 0) return iPoint;
- assert(iMc >= 0 && iMc < fpTofPoints->Size(iFile, iEvent));
+ assert(iMc >= 0 && iMc < fvpBrPoints[L1DetectorID::kTof]->Size(iFile, iEvent));
int iPointPrelim = fpMCData->FindInternalPointIndex(L1DetectorID::kTof, iMc, iEvent, iFile);
if (iPointPrelim == -1) { continue; }
iPoint = iPointPrelim;
@@ -369,54 +355,7 @@ int MCModule::MatchHitWithMc<L1DetectorID::kTof>(int iHitExt) const
return iPoint;
}
-// ---------------------------------------------------------------------------------------------------------------------
-//
-void MCModule::MatchPointsAndHits()
-{
- // Clear point indexes of each hit
- // NOTE: at the moment only the best point is stored to the hit, but we keep the possibility for keeping many
- // points for possible future module extension.
- for (auto& hit : (*fpvQaHits)) {
- hit.mcPointIds.clear();
- }
- // Match hit with point for different detectors
- auto DoMatch = [&](int iH, int iP) constexpr
- {
- if (iP > -1) {
- auto& hit = (*fpvQaHits)[iH];
- hit.mcPointIds.push_back_no_warning(iP);
- fpMCData->GetPoint(iP).AddHitID(iH);
- }
- };
-
- int iH = 0;
- if (fbUseMvd) {
- for (iH = (*fpvFstHitId)[int(L1DetectorID::kMvd)]; iH < (*fpvFstHitId)[int(L1DetectorID::kMvd) + 1]; ++iH) {
- DoMatch(iH, MatchHitWithMc<L1DetectorID::kMvd>((*fpvQaHits)[iH].ExtIndex));
- }
- }
- if (fbUseSts) {
- for (iH = (*fpvFstHitId)[int(L1DetectorID::kSts)]; iH < (*fpvFstHitId)[int(L1DetectorID::kSts) + 1]; ++iH) {
- DoMatch(iH, MatchHitWithMc<L1DetectorID::kSts>((*fpvQaHits)[iH].ExtIndex));
- }
- }
- if (fbUseSts) {
- for (iH = (*fpvFstHitId)[int(L1DetectorID::kMuch)]; iH < (*fpvFstHitId)[int(L1DetectorID::kMuch) + 1]; ++iH) {
- DoMatch(iH, MatchHitWithMc<L1DetectorID::kMuch>((*fpvQaHits)[iH].ExtIndex));
- }
- }
- if (fbUseTrd) {
- for (iH = (*fpvFstHitId)[int(L1DetectorID::kTrd)]; iH < (*fpvFstHitId)[int(L1DetectorID::kTrd) + 1]; ++iH) {
- DoMatch(iH, MatchHitWithMc<L1DetectorID::kTrd>((*fpvQaHits)[iH].ExtIndex));
- }
- }
- if (fbUseTof) {
- for (iH = (*fpvFstHitId)[int(L1DetectorID::kTof)]; iH < (*fpvFstHitId)[int(L1DetectorID::kTof) + 1]; ++iH) {
- DoMatch(iH, MatchHitWithMc<L1DetectorID::kTof>((*fpvQaHits)[iH].ExtIndex));
- }
- }
-}
// ---------------------------------------------------------------------------------------------------------------------
//
@@ -429,14 +368,15 @@ void MCModule::MatchRecoAndMCTracks()
auto& mNofHitsVsMCTrkID = trkRe.hitMap;
//mNofHitsVsMCTrkID.clear();
for (int iH : trkRe.Hits) {
- for (int iP : (*fpvQaHits)[iH].mcPointIds) {
- assert(iP > -1); // Should never be triggered
+ int iP = (*fpvQaHits)[iH].GetMCPointIndex();
+ L1_SHOW(iP);
+ if (iP > -1) {
int iTmc = fpMCData->GetPoint(iP).GetTrackId();
if (mNofHitsVsMCTrkID.find(iTmc) == mNofHitsVsMCTrkID.end()) { mNofHitsVsMCTrkID[iTmc] = 1; }
else {
mNofHitsVsMCTrkID[iTmc] += 1;
}
- } // loop over global point ids stored for hit: end
+ }
} // loop over hit ids stored for a reconstructed track: end
@@ -473,25 +413,6 @@ void MCModule::MatchRecoAndMCTracks()
} // loop over reconstructed tracks: end
}
-
-// ***********************
-// ** Accessors **
-// ***********************
-
-// ---------------------------------------------------------------------------------------------------------------------
-//
-void MCModule::SetDetector(L1DetectorID detID, bool flag)
-{
- switch (detID) {
- case L1DetectorID::kMvd: fbUseMvd = flag; break;
- case L1DetectorID::kSts: fbUseSts = flag; break;
- case L1DetectorID::kMuch: fbUseMuch = flag; break;
- case L1DetectorID::kTrd: fbUseTrd = flag; break;
- case L1DetectorID::kTof: fbUseTof = flag; break;
- case L1DetectorID::kEND: break;
- }
-}
-
// *******************************
// ** Utility functions **
// *******************************
@@ -521,43 +442,26 @@ void MCModule::CheckInit() const
if (!fpMCEventHeader) { throw std::logic_error("MC event header is unavailable"); }
// Check detectors initialization
- if (fbUseMvd) {
- if (!fpMvdPoints) { throw std::logic_error("MC points unavailable for MVD"); }
- if (!fpMvdHitMatches) { throw std::logic_error("Hit matches unavailable for MVD"); }
+ for (int iD = 0; iD < static_cast<int>(L1DetectorID::kEND); ++iD) {
+ if (fvbUseDet[iD]) {
+ if (!fvpBrPoints[iD]) { throw std::logic_error(Form("MC points are unavailable for %s", kDetName[iD])); }
+ if (!fvpBrHitMatches[iD]) { throw std::logic_error(Form("Hit matches are unavailable for %s", kDetName[iD])); }
+ }
}
- if (fbUseSts) {
- if (!fpStsPoints) { throw std::logic_error("MC points unavailable for STS"); }
- if (!fpStsHitMatches) { throw std::logic_error("Hit matches unavailable for STS"); }
+ if (fvbUseDet[L1DetectorID::kSts]) {
if (!fpStsClusterMatches) { throw std::logic_error("Cluster matches unavailable for STS"); }
if (!fpStsHits) { throw std::logic_error("Hits unavailable for STS"); }
}
-
- if (fbUseMuch) {
- if (!fpMuchPoints) { throw std::logic_error("MC points unavailable for MuCh"); }
- if (!fpMuchHitMatches) { throw std::logic_error("Hit matches unavailable for MuCh"); }
- }
-
- if (fbUseTrd) {
- if (!fpTrdPoints) { throw std::logic_error("MC points unavailable for TRD"); }
- if (!fpTrdHitMatches) { throw std::logic_error("Hit matches unavailable for TRD"); }
- }
-
- if (fbUseTof) {
- if (!fpTofPoints) { throw std::logic_error("MC points unavailable for TOF"); }
- if (!fpTofHitMatches) { throw std::logic_error("Hit matches unavailable for TOF"); }
- }
}
// ---------------------------------------------------------------------------------------------------------------------
//
template<>
-void MCModule::ReadMCPointsForDetector<L1DetectorID::kTof>(CbmMCDataArray* pPoints)
+void MCModule::ReadMCPointsForDetector<L1DetectorID::kTof>()
{
- for (const auto& key : fFileEventIDs) {
- int iFile = key.first;
- int iEvent = key.second;
-
+ auto* pPoints = fvpBrPoints[L1DetectorID::kTof];
+ for (const auto& [iFile, iEvent] : fFileEventIDs) {
// Prepare TODO
int nTofStationsGeo = fpParameters->GetNstationsGeometry(L1DetectorID::kTof);
std::vector<std::vector<int>> vTofPointToTrack(nTofStationsGeo); // TODO
@@ -573,7 +477,7 @@ void MCModule::ReadMCPointsForDetector<L1DetectorID::kTof>(CbmMCDataArray* pPoin
// Check whether a particular TOF point is matched to a hit
auto* pTofHits = dynamic_cast<TClonesArray*>(FairRootManager::Instance()->GetObject("TofHit"));
for (int iHit = 0; iHit < pTofHits->GetEntriesFast(); ++iHit) {
- auto* pHitMatch = L1_DYNAMIC_CAST<CbmMatch*>(fpTofHitMatches->At(iHit));
+ auto* pHitMatch = L1_DYNAMIC_CAST<CbmMatch*>(fvpBrHitMatches[L1DetectorID::kTof]->At(iHit));
LOG_IF(fatal, !pHitMatch) << "CA MC Module: hit match was not found for TOF hit " << iHit;
for (int iLink = 0; iLink < pHitMatch->GetNofLinks(); ++iLink) {
const auto& link = pHitMatch->GetLink(iLink);
@@ -627,14 +531,10 @@ void MCModule::ReadMCPointsForDetector<L1DetectorID::kTof>(CbmMCDataArray* pPoin
if (iStActive < 0) { continue; }
for (int iTrk = 0; iTrk < (int) vTofPointToTrack[iStLocGeo].size(); ++iTrk) {
if (vTofPointToTrack[iStLocGeo][iTrk] < 0) { continue; }
- ::ca::tools::MCPoint aPoint;
- if (FillMCPoint<L1DetectorID::kTof>(vTofPointToTrack[iStLocGeo][iTrk], iEvent, iFile, aPoint)) {
- aPoint.SetStationId(iStActive);
- int iExtGlob = fpMCData->GetPointGlobExtIndex(L1DetectorID::kTof, vTofPointToTrack[iStLocGeo][iTrk]);
- aPoint.SetExternalId(iExtGlob);
- int iPInt = fpMCData->GetNofPoints(); // assign an index of point in internal container
- if (aPoint.GetTrackId() > -1) { fpMCData->GetTrack(aPoint.GetTrackId()).AddPointIndex(iPInt); }
- fpMCData->AddPoint(aPoint);
+ auto oPoint = FillMCPoint<L1DetectorID::kTof>(vTofPointToTrack[iStLocGeo][iTrk], iEvent, iFile);
+ if (oPoint) {
+ oPoint->SetStationId(iStActive);
+ fpMCData->AddPoint(*oPoint);
}
} // loop over tracks: end
} // loop over geometry stations: end
@@ -648,34 +548,20 @@ void MCModule::ReadMCPoints()
int nPointsEstimated = 5 * fpMCData->GetNofTracks() * fpParameters->GetNstationsActive();
fpMCData->ReserveNofPoints(nPointsEstimated);
- int nPointsMvd = 0;
- int nPointsSts = 0;
- int nPointsMuch = 0;
- int nPointsTrd = 0;
- int nPointsTof = 0;
- for (const auto& key : fFileEventIDs) {
- int iFile = key.first;
- int iEvent = key.second;
-
- if (fbUseMvd) { nPointsMvd += fpMvdPoints->Size(iFile, iEvent); }
- if (fbUseSts) { nPointsSts += fpStsPoints->Size(iFile, iEvent); }
- if (fbUseMuch) { nPointsMuch += fpMuchPoints->Size(iFile, iEvent); }
- if (fbUseTrd) { nPointsTrd += fpTrdPoints->Size(iFile, iEvent); }
- if (fbUseTof) { nPointsTof += fpTofPoints->Size(iFile, iEvent); }
+ CbmCaDetIdArr_t<int> vNofPointsDet = {0};
+ for (const auto& [iFile, iEvent] : fFileEventIDs) {
+ for (int iD = 0; iD < static_cast<int>(vNofPointsDet.size()); ++iD) {
+ if (fvbUseDet[iD]) { vNofPointsDet[iD] = fvpBrPoints[iD]->Size(iFile, iEvent); }
+ fpMCData->SetNofPointsOrig(static_cast<L1DetectorID>(iD), vNofPointsDet[iD]);
+ }
}
- fpMCData->SetNofPointsOrig(L1DetectorID::kMvd, nPointsMvd);
- fpMCData->SetNofPointsOrig(L1DetectorID::kSts, nPointsSts);
- fpMCData->SetNofPointsOrig(L1DetectorID::kMuch, nPointsMuch);
- fpMCData->SetNofPointsOrig(L1DetectorID::kTrd, nPointsTrd);
- fpMCData->SetNofPointsOrig(L1DetectorID::kTof, nPointsTof);
-
// ----- Read MC points in MVD, STS, MuCh, TRD and TOF
- if (fbUseMvd) { this->ReadMCPointsForDetector<L1DetectorID::kMvd>(fpMvdPoints); }
- if (fbUseSts) { this->ReadMCPointsForDetector<L1DetectorID::kSts>(fpStsPoints); }
- if (fbUseMuch) { this->ReadMCPointsForDetector<L1DetectorID::kMuch>(fpMuchPoints); }
- if (fbUseTrd) { this->ReadMCPointsForDetector<L1DetectorID::kTrd>(fpTrdPoints); }
- if (fbUseTof) { this->ReadMCPointsForDetector<L1DetectorID::kTof>(fpTofPoints); }
+ if (fvbUseDet[L1DetectorID::kMvd]) { this->ReadMCPointsForDetector<L1DetectorID::kMvd>(); }
+ if (fvbUseDet[L1DetectorID::kSts]) { this->ReadMCPointsForDetector<L1DetectorID::kSts>(); }
+ if (fvbUseDet[L1DetectorID::kMuch]) { this->ReadMCPointsForDetector<L1DetectorID::kMuch>(); }
+ if (fvbUseDet[L1DetectorID::kTrd]) { this->ReadMCPointsForDetector<L1DetectorID::kTrd>(); }
+ if (fvbUseDet[L1DetectorID::kTof]) { this->ReadMCPointsForDetector<L1DetectorID::kTof>(); }
}
// ---------------------------------------------------------------------------------------------------------------------
diff --git a/reco/L1/CbmCaMCModule.h b/reco/L1/CbmCaMCModule.h
index 95e5ecf16735fb19389169469e2d6fdc81948352..e1f2b90296fc06f8cb9c8f782f2717e008eb2536 100644
--- a/reco/L1/CbmCaMCModule.h
+++ b/reco/L1/CbmCaMCModule.h
@@ -29,6 +29,7 @@
#include "TDatabasePDG.h"
#include <numeric>
+#include <optional>
#include <string>
#include <string_view>
#include <type_traits>
@@ -52,6 +53,8 @@ namespace cbm::ca
/// @brief Class CbmCaPerformance is an interface to communicate between
///
class MCModule {
+
+
public:
/// @brief Constructor
/// @param verbosity Verbosity level
@@ -73,17 +76,13 @@ namespace cbm::ca
/// @brief Move assignment operator
MCModule& operator=(MCModule&&) = delete;
- /// @brief Creates hits from MC points
- ///
- /// This function creates hits from MC points in selected tracking stations and re-fills hit arrays.
- void CreateHitsFromPoints() {}
-
- /// @brief Sets hit parameters from the matched MC point
- void AdjustHitsToPoints();
/// @brief Defines performance action in the end of the run
void Finish();
+ /// @brief Gets a pointer to MC data object
+ const ::ca::tools::MCData* GetMCData() const { return fpMCData; }
+
/// @brief Defines performance action in the beginning of each event or time slice
/// @note This function should be called before hits initialization
/// @param pEvent Pointer to a current CbmEvent
@@ -93,6 +92,16 @@ namespace cbm::ca
/// @return Success flag
bool InitRun();
+ /// @brief Matches hit with MC point
+ /// @tparam DetId Detector ID
+ /// @param iHitExt External index of hit
+ /// @return MC-point index in fvMCPoints array
+ ///
+ /// This method finds a match for a given hit or matches for hits clusters (in case of STS), finds the best
+ /// link in the match and returns the corresponding global ID of the MC points.
+ template<L1DetectorID DetId>
+ int MatchHitWithMc(int iHitExt) const;
+
/// @brief Match reconstructed and MC data
///
/// Runs matching of hits with points and reconstructed tracks with MC ones. Reconstructed tracks are updated with
@@ -115,7 +124,7 @@ namespace cbm::ca
/// @param detID Id of detector
/// @param flag Flag: true - detector is used
/// @note Should be called before this->Init()
- void SetDetector(L1DetectorID detID, bool flag);
+ void SetDetector(L1DetectorID detID, bool flag) { fvbUseDet[detID] = flag; }
/// @brief Sets legacy event mode:
/// @param flag Flag:
@@ -123,26 +132,6 @@ namespace cbm::ca
/// - false: runs on time-slice base
void SetLegacyEventMode(bool flag) { fbLegacyEventMode = flag; }
- /// @brief Creates hits from MC points for a particular detector and station
- /// @param detID DetectorID
- /// @param iStLoc Local index of station, provided by geometry
- /// @param du Error u-coordinate [cm]
- /// @param dv Error v-coordinate [cm]
- /// @param dt Error of time [ns]
- /// @param ifSmear Flag:
- /// - true: MC point quantities are smeared along the error by gaussian,
- /// - false: Precise point quantities are used
- ///void SetHitsFromPoints(L1DetectorID detID, int iStLoc, double du, double dv, double dt, bool ifSmear) {}
-
- /// @brief Adjusts existing reconstructed hits to MC points for a particular detector and station
- /// @param detID DetectorID
- /// @param iStLoc Local index of station, provided by geometry
- /// @param ifSmear Flag:
- /// - true: MC point quantities are smeared along the error by gaussian,
- /// - false: Precise point quantities are used
- ///void SetHitsFromPoints(L1DetectorID detID, int iStLoc, bool ifSmear) {}
-
-
/// @brief Registers MC data object
/// @param mcData Instance of MC data
void RegisterMCData(::ca::tools::MCData& mcData) { fpMCData = &mcData; }
@@ -175,20 +164,11 @@ namespace cbm::ca
/// consecutive stations containing a hit or point and number of stations and points with hits.
void InitTrackInfo();
-
- /// @brief Matches hit with MC point
- /// @tparam DetId Detector ID
- /// @param iHitExt External index of hit
- /// @return MC-point index in fvMCPoints array
- ///
- /// This method finds a match for a given hit or matches for hits clusters (in case of STS), finds the best
- /// link in the match and returns the corresponding global ID of the MC points.
- template<L1DetectorID DetId>
- int MatchHitWithMc(int iHitExt) const;
-
- /// @brief Match MC points and reconstructed hits
+ /// @brief Match sets of MC points and reconstructed hits for a given detector
+ /// @tparam DetID Index of the detector
///
/// Writes indexes of MC points to each hit and indexes of hits to each MC point.
+ template<L1DetectorID DetID>
void MatchPointsAndHits();
/// @brief Matches reconstructed tracks with MC tracks
@@ -209,7 +189,7 @@ namespace cbm::ca
/// @brief Reads MC points in particular detector
template<L1DetectorID DetID>
- void ReadMCPointsForDetector(CbmMCDataArray* pPoints);
+ void ReadMCPointsForDetector();
/// @brief Fills a single detector-specific MC point
/// @tparam DetID Detector subsystem ID
@@ -220,42 +200,23 @@ namespace cbm::ca
/// @return true Point is correct and is to be saved to the MCData object
/// @return false Point is incorrect and will be ignored
template<L1DetectorID DetID>
- bool FillMCPoint(int iExtId, int iEvent, int iFile, ::ca::tools::MCPoint& point);
-
- std::shared_ptr<L1Parameters> fpParameters = nullptr; ///< Pointer to tracking parameters object
+ std::optional<::ca::tools::MCPoint> FillMCPoint(int iExtId, int iEvent, int iFile);
// ------ Flags
- bool fbUseMvd = false; ///< is MVD used
- bool fbUseSts = false; ///< is STS used
- bool fbUseMuch = false; ///< is MuCh used
- bool fbUseTrd = false; ///< is TRD used
- bool fbUseTof = false; ///< is TOF used
-
- bool fbLegacyEventMode = false; ///< if tracking uses events instead of time-slices (back compatibility)
- int fVerbose = 0; ///< Verbosity level
- int fPerformanceMode = -1; ///< Mode of performance
+ CbmCaDetIdArr_t<bool> fvbUseDet = {{false}}; ///< Flag: is detector subsystem used
+ bool fbLegacyEventMode = false; ///< if tracking uses events instead of time-slices (back compatibility)
+ int fVerbose = 0; ///< Verbosity level
+ int fPerformanceMode = -1; ///< Mode of performance
+ std::shared_ptr<L1Parameters> fpParameters = nullptr; ///< Pointer to tracking parameters object
// ------ Input data branches
- const CbmTimeSlice* fpTimeSlice = nullptr; ///< Current time slice
-
- // Mc-event
+ const CbmTimeSlice* fpTimeSlice = nullptr; ///< Current time slice
CbmMCEventList* fpMCEventList = nullptr; ///< MC event list
CbmMCDataObject* fpMCEventHeader = nullptr; ///< MC event header
CbmMCDataArray* fpMCTracks = nullptr; ///< MC tracks input
- // Mc-points
- CbmMCDataArray* fpMvdPoints = nullptr; ///< MVD MC-points input container
- CbmMCDataArray* fpStsPoints = nullptr; ///< STS MC-points input container
- CbmMCDataArray* fpMuchPoints = nullptr; ///< MuCh MC-points input container
- CbmMCDataArray* fpTrdPoints = nullptr; ///< TRD MC-points input container
- CbmMCDataArray* fpTofPoints = nullptr; ///< TOF MC-points input container
-
- // Matches
- TClonesArray* fpMvdHitMatches = nullptr; ///< Array of MVD hit matches
- TClonesArray* fpStsHitMatches = nullptr; ///< Array of STS hit matches
- TClonesArray* fpMuchHitMatches = nullptr; ///< Array of MuCh hit matches
- TClonesArray* fpTrdHitMatches = nullptr; ///< Array of TRD hit matches
- TClonesArray* fpTofHitMatches = nullptr; ///< Array of TOF hit matches
+ CbmCaDetIdArr_t<CbmMCDataArray*> fvpBrPoints = {{nullptr}}; ///< Array of points vs. detector
+ CbmCaDetIdArr_t<TClonesArray*> fvpBrHitMatches = {{nullptr}}; ///< Array of hit match branches vs. detector
TClonesArray* fpStsHits = nullptr; ///< Array of STS hits (currently needed for matching)
TClonesArray* fpStsClusterMatches = nullptr; ///< Array of STS cluster matches
@@ -263,14 +224,15 @@ namespace cbm::ca
// Matching information
std::set<std::pair<int, int>> fFileEventIDs; ///< Set of file and event indexes: first - iFile, second - iEvent
- // MC data
+ // ----- Internal MC data
::ca::tools::MCData* fpMCData = nullptr; ///< MC information (hits and tracks) instance
- // Reconstructed data
+ // ----- Internal reconstructed data
L1Vector<CbmL1Track>* fpvRecoTracks = nullptr; ///< Pointer to reconstructed track container
L1Vector<CbmL1HitId>* fpvHitIds = nullptr; ///< Pointer to hit index container
L1Vector<CbmL1HitDebugInfo>* fpvQaHits = nullptr; ///< Pointer to QA hit container
+
/// @brief Pointer to array of first hit indexes in the detector subsystem
///
/// This array must be initialized in the run initialization function.
@@ -286,8 +248,10 @@ namespace cbm::ca
// ---------------------------------------------------------------------------------------------------------------------
//
template<L1DetectorID DetID>
-bool cbm::ca::MCModule::FillMCPoint(int iExtId, int iEvent, int iFile, ::ca::tools::MCPoint& point)
+std::optional<::ca::tools::MCPoint> cbm::ca::MCModule::FillMCPoint(int iExtId, int iEvent, int iFile)
{
+ auto oPoint = std::make_optional<::ca::tools::MCPoint>();
+
// ----- Get positions, momenta, time and track ID
TVector3 posIn; // Position at entrance to station [cm]
TVector3 posOut; // Position at exist of station [cm]
@@ -295,13 +259,15 @@ bool cbm::ca::MCModule::FillMCPoint(int iExtId, int iEvent, int iFile, ::ca::too
TVector3 momOut; // 3-momentum at exit of station [GeV/c]
double time = undef::kD; // Time of MC point (after beginning of the event) [ns]
int iTmcExt = undef::kI32; // Track ID in the external container
- // MVD
+
+ auto* pBrPoints = fvpBrPoints[DetID];
+
if constexpr (L1DetectorID::kMvd == DetID) {
- auto* pExtPoint = dynamic_cast<CbmMvdPoint*>(fpMvdPoints->Get(iFile, iEvent, iExtId));
+ auto* pExtPoint = dynamic_cast<CbmMvdPoint*>(pBrPoints->Get(iFile, iEvent, iExtId));
if (!pExtPoint) {
LOG(warn) << "CbmCaMCModule: MVD MC point with iExtId = " << iExtId << ", iEvent = " << iEvent
<< ", iFile = " << iFile << " does not exist";
- return false;
+ return std::nullopt;
}
pExtPoint->Position(posIn);
pExtPoint->PositionOut(posOut);
@@ -312,11 +278,11 @@ bool cbm::ca::MCModule::FillMCPoint(int iExtId, int iEvent, int iFile, ::ca::too
}
// STS
else if constexpr (L1DetectorID::kSts == DetID) {
- auto* pExtPoint = dynamic_cast<CbmStsPoint*>(fpStsPoints->Get(iFile, iEvent, iExtId));
+ auto* pExtPoint = dynamic_cast<CbmStsPoint*>(pBrPoints->Get(iFile, iEvent, iExtId));
if (!pExtPoint) {
LOG(warn) << "CbmCaMCModule: STS MC point with iExtId = " << iExtId << ", iEvent = " << iEvent
<< ", iFile = " << iFile << " does not exist";
- return false;
+ return std::nullopt;
}
pExtPoint->Position(posIn);
pExtPoint->PositionOut(posOut);
@@ -327,11 +293,11 @@ bool cbm::ca::MCModule::FillMCPoint(int iExtId, int iEvent, int iFile, ::ca::too
}
// MuCh
else if constexpr (L1DetectorID::kMuch == DetID) {
- auto* pExtPoint = dynamic_cast<CbmMuchPoint*>(fpMuchPoints->Get(iFile, iEvent, iExtId));
+ auto* pExtPoint = dynamic_cast<CbmMuchPoint*>(pBrPoints->Get(iFile, iEvent, iExtId));
if (!pExtPoint) {
LOG(warn) << "CbmCaMCModule: MuCh MC point with iExtId = " << iExtId << ", iEvent = " << iEvent
<< ", iFile = " << iFile << " does not exist";
- return false;
+ return std::nullopt;
}
pExtPoint->Position(posIn);
pExtPoint->PositionOut(posOut);
@@ -342,11 +308,11 @@ bool cbm::ca::MCModule::FillMCPoint(int iExtId, int iEvent, int iFile, ::ca::too
}
// TRD
else if constexpr (L1DetectorID::kTrd == DetID) {
- auto* pExtPoint = dynamic_cast<CbmTrdPoint*>(fpTrdPoints->Get(iFile, iEvent, iExtId));
+ auto* pExtPoint = dynamic_cast<CbmTrdPoint*>(pBrPoints->Get(iFile, iEvent, iExtId));
if (!pExtPoint) {
LOG(warn) << "CbmCaMCModule: TRD MC point with iExtId = " << iExtId << ", iEvent = " << iEvent
<< ", iFile = " << iFile << " does not exist";
- return false;
+ return std::nullopt;
}
pExtPoint->Position(posIn);
pExtPoint->PositionOut(posOut);
@@ -357,11 +323,11 @@ bool cbm::ca::MCModule::FillMCPoint(int iExtId, int iEvent, int iFile, ::ca::too
}
// TOF
else if constexpr (L1DetectorID::kTof == DetID) {
- auto* pExtPoint = dynamic_cast<CbmTofPoint*>(fpTofPoints->Get(iFile, iEvent, iExtId));
+ auto* pExtPoint = dynamic_cast<CbmTofPoint*>(pBrPoints->Get(iFile, iEvent, iExtId));
if (!pExtPoint) {
LOG(warn) << "CbmCaMCModule: TOF MC point with iExtId = " << iExtId << ", iEvent = " << iEvent
<< ", iFile = " << iFile << " does not exist";
- return false;
+ return std::nullopt;
}
pExtPoint->Position(posIn);
pExtPoint->Position(posOut);
@@ -373,7 +339,7 @@ bool cbm::ca::MCModule::FillMCPoint(int iExtId, int iEvent, int iFile, ::ca::too
if (iTmcExt < 0) {
LOG(warn) << "CbmCaMCModule: For MC point with iExtId = " << iExtId << ", iEvent = " << iEvent
<< ", iFile = " << iFile << " MC track is undefined (has ID = " << iTmcExt << ')';
- return false;
+ return std::nullopt;
}
TVector3 posMid = 0.5 * (posIn + posOut);
TVector3 momMid = 0.5 * (momIn + momOut);
@@ -421,79 +387,91 @@ bool cbm::ca::MCModule::FillMCPoint(int iExtId, int iEvent, int iFile, ::ca::too
LOG(warn) << "CbmCaMCModule: MC point with iExtId = " << iExtId << ", iEvent = " << iEvent
<< ", iFile = " << iFile << " and det id " << int(DetID) << " fell out of the TS duration [" << startT
<< ", " << endT << "] with measured time = " << time << " [ns]";
- return false;
+ return std::nullopt;
}
}
// ----- Fill MC point
- point.SetExternalId(iExtId);
- point.SetEventId(iEvent);
- point.SetFileId(iFile);
- point.SetTime(time);
- point.SetX(posMid.X());
- point.SetY(posMid.Y());
- point.SetZ(posMid.Z());
- point.SetXIn(posIn.X());
- point.SetYIn(posIn.Y());
- point.SetZIn(posIn.Z());
- point.SetXOut(posOut.X());
- point.SetYOut(posOut.Y());
- point.SetZOut(posOut.Z());
- point.SetPx(momMid.X());
- point.SetPy(momMid.Y());
- point.SetPz(momMid.Z());
- point.SetPxIn(momIn.X());
- point.SetPyIn(momIn.Y());
- point.SetPzIn(momIn.Z());
- point.SetPxOut(momOut.X());
- point.SetPyOut(momOut.Y());
- point.SetPzOut(momOut.Z());
-
+ oPoint->SetExternalId(fpMCData->GetPointGlobExtIndex(DetID, iExtId));
+ oPoint->SetEventId(iEvent);
+ oPoint->SetFileId(iFile);
+ oPoint->SetTime(time);
+ oPoint->SetX(posMid.X());
+ oPoint->SetY(posMid.Y());
+ oPoint->SetZ(posMid.Z());
+ oPoint->SetXIn(posIn.X());
+ oPoint->SetYIn(posIn.Y());
+ oPoint->SetZIn(posIn.Z());
+ oPoint->SetXOut(posOut.X());
+ oPoint->SetYOut(posOut.Y());
+ oPoint->SetZOut(posOut.Z());
+ oPoint->SetPx(momMid.X());
+ oPoint->SetPy(momMid.Y());
+ oPoint->SetPz(momMid.Z());
+ oPoint->SetPxIn(momIn.X());
+ oPoint->SetPyIn(momIn.Y());
+ oPoint->SetPzIn(momIn.Z());
+ oPoint->SetPxOut(momOut.X());
+ oPoint->SetPyOut(momOut.Y());
+ oPoint->SetPzOut(momOut.Z());
+
+ // Select current number of points as a local id of points
+ oPoint->SetId(fpMCData->GetNofPoints());
+
+ // Match MC track and point to each other
int iTmcInt = fpMCData->FindInternalTrackIndex(iTmcExt, iEvent, iFile);
- point.SetId(fpMCData->GetNofPoints()); // select current number of points as a local id of points
- point.SetTrackId(iTmcInt);
- point.SetStationId(stationID);
- point.SetDetectorId(DetID);
+ oPoint->SetTrackId(iTmcInt);
+ if (iTmcInt > -1) { fpMCData->GetTrack(iTmcInt).AddPointIndex(oPoint->GetId()); }
+
+ oPoint->SetStationId(stationID);
+ oPoint->SetDetectorId(DetID);
auto* pExtTrk = L1_DYNAMIC_CAST<CbmMCTrack*>(fpMCTracks->Get(iFile, iEvent, iTmcExt));
if (!pExtTrk) {
LOG(warn) << "CbmCaMCModule: MC track with iTmcExt = " << iTmcExt << ", iEvent = " << iEvent
<< ", iFile = " << iFile << " MC track is undefined (nullptr)";
- return false;
+ return std::nullopt;
}
- point.SetPdgCode(pExtTrk->GetPdgCode());
- point.SetMotherId(pExtTrk->GetMotherId());
+ oPoint->SetPdgCode(pExtTrk->GetPdgCode());
+ oPoint->SetMotherId(pExtTrk->GetMotherId());
- auto* pPdgDB = TDatabasePDG::Instance()->GetParticle(point.GetPdgCode());
- point.SetMass(pPdgDB ? pPdgDB->Mass() : 0.); /// TODO: Set from track
- point.SetCharge(pPdgDB ? pPdgDB->Charge() / 3. : 0.);
+ auto* pPdgDB = TDatabasePDG::Instance()->GetParticle(oPoint->GetPdgCode());
+ oPoint->SetMass(pPdgDB ? pPdgDB->Mass() : 0.); /// TODO: Set from track
+ oPoint->SetCharge(pPdgDB ? pPdgDB->Charge() / 3. : 0.);
- return true;
+ return oPoint;
}
// ---------------------------------------------------------------------------------------------------------------------
-// NOTE: template is used, because another template function FillMCPoint is used inside
+//
template<L1DetectorID DetID>
-void cbm::ca::MCModule::ReadMCPointsForDetector(CbmMCDataArray* pPoints)
+void cbm::ca::MCModule::MatchPointsAndHits()
{
- if constexpr (L1DetectorID::kTof != DetID) {
- for (const auto& key : fFileEventIDs) {
- int iFile = key.first;
- int iEvent = key.second;
- int nPointsEvent = pPoints->Size(iFile, iEvent);
- for (int iP = 0; iP < nPointsEvent; ++iP) {
- ::ca::tools::MCPoint aPoint;
- if (FillMCPoint<DetID>(iP, iEvent, iFile, aPoint)) {
- aPoint.SetExternalId(fpMCData->GetPointGlobExtIndex(DetID, iP));
- int iPInt = fpMCData->GetNofPoints(); // assign an index of point in internal container
- if (aPoint.GetTrackId() > -1) { fpMCData->GetTrack(aPoint.GetTrackId()).AddPointIndex(iPInt); }
- fpMCData->AddPoint(aPoint);
- }
- } // iP: end
- } // key: end
+ if (!fvbUseDet[DetID]) { return; }
+
+ for (int iH = (*fpvFstHitId)[static_cast<int>(DetID)]; iH < (*fpvFstHitId)[static_cast<int>(DetID) + 1]; ++iH) {
+ auto& hit = (*fpvQaHits)[iH];
+ int iP = MatchHitWithMc<DetID>(hit.ExtIndex);
+ if (iP > -1) {
+ hit.SetMCPointIndex(iP);
+ fpMCData->GetPoint(iP).AddHitID(iH);
+ }
}
}
+// ---------------------------------------------------------------------------------------------------------------------
+// NOTE: template is used, because another template function FillMCPoint is used inside
+template<L1DetectorID DetID>
+void cbm::ca::MCModule::ReadMCPointsForDetector()
+{
+ for (const auto& [iFile, iEvent] : fFileEventIDs) {
+ int nPointsEvent = fvpBrPoints[DetID]->Size(iFile, iEvent);
+ for (int iP = 0; iP < nPointsEvent; ++iP) {
+ std::optional<::ca::tools::MCPoint> oPoint = FillMCPoint<DetID>(iP, iEvent, iFile);
+ if (oPoint) { fpMCData->AddPoint(*oPoint); }
+ } // iP: end
+ } // key: end
+}
#endif // CbmCaMCModule_h
diff --git a/reco/L1/CbmCaMCModule.h.orig b/reco/L1/CbmCaMCModule.h.orig
new file mode 100644
index 0000000000000000000000000000000000000000..ebe5bbe25f7c0288b8c5f4ce10ffefeadc2a4266
--- /dev/null
+++ b/reco/L1/CbmCaMCModule.h.orig
@@ -0,0 +1,503 @@
+/* Copyright (C) 2022-2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// @file CbmCaMCModule.h
+/// @brief CA Tracking performance interface for CBM (header)
+/// @since 23.09.2022
+/// @author S.Zharko <s.zharko@gsi.de>
+
+#ifndef CbmCaMCModule_h
+#define CbmCaMCModule_h 1
+
+#include "CbmL1DetectorID.h"
+#include "CbmL1Hit.h"
+#include "CbmL1Track.h"
+#include "CbmLink.h"
+#include "CbmMCDataArray.h"
+#include "CbmMCEventList.h"
+#include "CbmMCTrack.h"
+#include "CbmMatch.h"
+#include "CbmMuchPoint.h"
+#include "CbmMvdPoint.h"
+#include "CbmStsPoint.h"
+#include "CbmTimeSlice.h"
+#include "CbmTofPoint.h"
+#include "CbmTrdPoint.h"
+
+#include "TClonesArray.h"
+#include "TDatabasePDG.h"
+
+#include <numeric>
+#include <string>
+#include <string_view>
+#include <type_traits>
+
+#include "CaToolsMCData.h"
+#include "CaToolsMCPoint.h"
+#include "L1Constants.h"
+#include "L1InputData.h"
+#include "L1Parameters.h"
+#include "L1Undef.h"
+
+class CbmEvent;
+class CbmMCDataObject;
+class CbmL1HitDebugInfo;
+class CbmL1Track;
+
+enum class L1DetectorID;
+
+namespace cbm::ca
+{
+<<<<<<< HEAD
+ /// @brief Class CbmCaPerformance is an interface to communicate between
+=======
+ /// Class CbmCaPerformance is an interface to communicate between
+>>>>>>> CA: QA updates
+ ///
+ class MCModule {
+ public:
+ /// @brief Constructor
+ /// @param verbosity Verbosity level
+ /// @param perfMode Performance mode (defines cut on number of consecutive stations with hit or point)
+ MCModule(int verb = 0, int perfMode = 3) : fVerbose(verb), fPerformanceMode(perfMode) {}
+
+ /// @brief Destructor
+ ~MCModule() = default;
+
+ /// @brief Copy constructor
+ MCModule(const MCModule&) = delete;
+
+ /// @brief Move constructor
+ MCModule(MCModule&&) = delete;
+
+ /// @brief Copy assignment operator
+ MCModule& operator=(const MCModule&) = delete;
+
+ /// @brief Move assignment operator
+ MCModule& operator=(MCModule&&) = delete;
+
+ /// @brief Creates hits from MC points
+ ///
+ /// This function creates hits from MC points in selected tracking stations and re-fills hit arrays.
+ void CreateHitsFromPoints() {}
+
+ /// @brief Sets hit parameters from the matched MC point
+ void AdjustHitsToPoints();
+
+ /// @brief Defines performance action in the end of the run
+ void Finish();
+
+ /// @brief Defines performance action in the beginning of each event or time slice
+ /// @note This function should be called before hits initialization
+ /// @param pEvent Pointer to a current CbmEvent
+ void InitEvent(CbmEvent* pEvent);
+
+ /// @brief Defines action on the module in the beginning of the run
+ /// @return Success flag
+ bool InitRun();
+
+ /// @brief Match reconstructed and MC data
+ ///
+ /// Runs matching of hits with points and reconstructed tracks with MC ones. Reconstructed tracks are updated with
+ /// true information.
+ void MatchRecoAndMC();
+
+ /// @brief Processes event
+ ///
+ /// Fills histograms and tables, should be called after the tracking done
+ void ProcessEvent(CbmEvent* pEvent);
+
+ /// @brief Sets first hit indexes container in different detectors
+ /// @param source Array of indexes
+ void RegisterFirstHitIndexes(const std::array<int, L1Constants::size::kMaxNdetectors + 1>& source)
+ {
+ fpvFstHitId = &source;
+ }
+
+ /// @brief Sets used detector subsystems
+ /// @param detID Id of detector
+ /// @param flag Flag: true - detector is used
+ /// @note Should be called before this->Init()
+ void SetDetector(L1DetectorID detID, bool flag);
+
+ /// @brief Sets legacy event mode:
+ /// @param flag Flag:
+ /// - true: runs on events base,
+ /// - false: runs on time-slice base
+ void SetLegacyEventMode(bool flag) { fbLegacyEventMode = flag; }
+
+ /// @brief Creates hits from MC points for a particular detector and station
+ /// @param detID DetectorID
+ /// @param iStLoc Local index of station, provided by geometry
+ /// @param du Error u-coordinate [cm]
+ /// @param dv Error v-coordinate [cm]
+ /// @param dt Error of time [ns]
+ /// @param ifSmear Flag:
+ /// - true: MC point quantities are smeared along the error by gaussian,
+ /// - false: Precise point quantities are used
+ ///void SetHitsFromPoints(L1DetectorID detID, int iStLoc, double du, double dv, double dt, bool ifSmear) {}
+
+ /// @brief Adjusts existing reconstructed hits to MC points for a particular detector and station
+ /// @param detID DetectorID
+ /// @param iStLoc Local index of station, provided by geometry
+ /// @param ifSmear Flag:
+ /// - true: MC point quantities are smeared along the error by gaussian,
+ /// - false: Precise point quantities are used
+ ///void SetHitsFromPoints(L1DetectorID detID, int iStLoc, bool ifSmear) {}
+
+
+ /// @brief Registers MC data object
+ /// @param mcData Instance of MC data
+ void RegisterMCData(::ca::tools::MCData& mcData) { fpMCData = &mcData; }
+
+ /// @brief Registers reconstructed track container
+ /// @param vRecoTrack Reference to reconstructed track container
+ void RegisterRecoTrackContainer(L1Vector<CbmL1Track>& vRecoTracks) { fpvRecoTracks = &vRecoTracks; }
+
+ /// @brief Registers hit index container
+ /// @param vHitIds Reference to hit index container
+ void RegisterHitIndexContainer(L1Vector<CbmL1HitId>& vHitIds) { fpvHitIds = &vHitIds; }
+
+ /// @brief Registers CA parameters object
+ /// @param pParameters A shared pointer to the parameters object
+ void RegisterParameters(std::shared_ptr<L1Parameters>& pParameters) { fpParameters = pParameters; }
+
+ /// @brief Registers debug hit container
+ /// @param vQaHits Reference to debug hit container
+ void RegisterQaHitContainer(L1Vector<CbmL1HitDebugInfo>& vQaHits) { fpvQaHits = &vQaHits; }
+
+ private:
+ /// @brief Check class initialization
+ /// @note The function throws std::logic_error, if initialization is incomplete
+ void CheckInit() const;
+
+ /// @brief Initializes MC track
+ ///
+ /// Initializes information about arrangement of points and hits of MC tracks within stations and the status
+ /// of track ability to be reconstructed, calculates max number of points and hits on a station, number of
+ /// consecutive stations containing a hit or point and number of stations and points with hits.
+ void InitTrackInfo();
+
+
+ /// @brief Matches hit with MC point
+ /// @tparam DetId Detector ID
+ /// @param iHitExt External index of hit
+ /// @return MC-point index in fvMCPoints array
+ ///
+ /// This method finds a match for a given hit or matches for hits clusters (in case of STS), finds the best
+ /// link in the match and returns the corresponding global ID of the MC points.
+ template<L1DetectorID DetId>
+ int MatchHitWithMc(int iHitExt) const;
+
+ /// @brief Match MC points and reconstructed hits
+ ///
+ /// Writes indexes of MC points to each hit and indexes of hits to each MC point.
+ void MatchPointsAndHits();
+
+ /// @brief Matches reconstructed tracks with MC tracks
+ ///
+ /// In the procedure, the maps of associated MC track indexes to corresponding number of hits are filled out and the
+ /// max purity is calculated for each reconstructed track in the TS/event. If the value of purity for a given MC track
+ /// is larger then a threshold, the corresponding track index is saved to the reconstructed track object, in the same
+ /// time the index of the reconstructed track object is saved to this MC track object. If purity for the MC track does
+ /// not pass the threshold, its index will not be accounted as a matched track, and the index of reconstructed track
+ /// will be added as an index of "touched" track.
+ void MatchRecoAndMCTracks();
+
+ /// @brief Reads MC tracks from external trees and saves them to MCDataObject
+ void ReadMCTracks();
+
+ /// @brief Reads MC points from external trees and saves them to MCDataObject
+ void ReadMCPoints();
+
+ /// @brief Reads MC points in particular detector
+ template<L1DetectorID DetID>
+ void ReadMCPointsForDetector(CbmMCDataArray* pPoints);
+
+ /// @brief Fills a single detector-specific MC point
+ /// @tparam DetID Detector subsystem ID
+ /// @param[in] iExtId Index of point in the external points container
+ /// @param[in] iEvent EventID of point in the external points container
+ /// @param[in] iFile FileID of point int the external points container
+ /// @param[out] intMCPoint Reference to the internal tracking MC point object (ca::tools::MCData)
+ /// @return true Point is correct and is to be saved to the MCData object
+ /// @return false Point is incorrect and will be ignored
+ template<L1DetectorID DetID>
+ bool FillMCPoint(int iExtId, int iEvent, int iFile, ::ca::tools::MCPoint& point);
+
+ std::shared_ptr<L1Parameters> fpParameters = nullptr; ///< Pointer to tracking parameters object
+
+ // ------ Flags
+ bool fbUseMvd = false; ///< is MVD used
+ bool fbUseSts = false; ///< is STS used
+ bool fbUseMuch = false; ///< is MuCh used
+ bool fbUseTrd = false; ///< is TRD used
+ bool fbUseTof = false; ///< is TOF used
+
+ bool fbLegacyEventMode = false; ///< if tracking uses events instead of time-slices (back compatibility)
+ int fVerbose = 0; ///< Verbosity level
+ int fPerformanceMode = -1; ///< Mode of performance
+
+ // ------ Input data branches
+ const CbmTimeSlice* fpTimeSlice = nullptr; ///< Current time slice
+
+ // Mc-event
+ CbmMCEventList* fpMCEventList = nullptr; ///< MC event list
+ CbmMCDataObject* fpMCEventHeader = nullptr; ///< MC event header
+ CbmMCDataArray* fpMCTracks = nullptr; ///< MC tracks input
+
+ // Mc-points
+ CbmMCDataArray* fpMvdPoints = nullptr; ///< MVD MC-points input container
+ CbmMCDataArray* fpStsPoints = nullptr; ///< STS MC-points input container
+ CbmMCDataArray* fpMuchPoints = nullptr; ///< MuCh MC-points input container
+ CbmMCDataArray* fpTrdPoints = nullptr; ///< TRD MC-points input container
+ CbmMCDataArray* fpTofPoints = nullptr; ///< TOF MC-points input container
+
+ // Matches
+ TClonesArray* fpMvdHitMatches = nullptr; ///< Array of MVD hit matches
+ TClonesArray* fpStsHitMatches = nullptr; ///< Array of STS hit matches
+ TClonesArray* fpMuchHitMatches = nullptr; ///< Array of MuCh hit matches
+ TClonesArray* fpTrdHitMatches = nullptr; ///< Array of TRD hit matches
+ TClonesArray* fpTofHitMatches = nullptr; ///< Array of TOF hit matches
+
+ TClonesArray* fpStsHits = nullptr; ///< Array of STS hits (currently needed for matching)
+ TClonesArray* fpStsClusterMatches = nullptr; ///< Array of STS cluster matches
+
+ // Matching information
+ std::set<std::pair<int, int>> fFileEventIDs; ///< Set of file and event indexes: first - iFile, second - iEvent
+
+ // MC data
+ ::ca::tools::MCData* fpMCData = nullptr; ///< MC information (hits and tracks) instance
+
+ // Reconstructed data
+ L1Vector<CbmL1Track>* fpvRecoTracks = nullptr; ///< Pointer to reconstructed track container
+ L1Vector<CbmL1HitId>* fpvHitIds = nullptr; ///< Pointer to hit index container
+ L1Vector<CbmL1HitDebugInfo>* fpvQaHits = nullptr; ///< Pointer to QA hit container
+
+ /// @brief Pointer to array of first hit indexes in the detector subsystem
+ ///
+ /// This array must be initialized in the run initialization function.
+ const std::array<int, L1Constants::size::kMaxNdetectors + 1>* fpvFstHitId = nullptr;
+ };
+} // namespace cbm::ca
+
+
+// **********************************************
+// ** Template function implementation **
+// **********************************************
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+template<L1DetectorID DetID>
+bool cbm::ca::MCModule::FillMCPoint(int iExtId, int iEvent, int iFile, ::ca::tools::MCPoint& point)
+{
+ // ----- Get positions, momenta, time and track ID
+ TVector3 posIn; // Position at entrance to station [cm]
+ TVector3 posOut; // Position at exist of station [cm]
+ TVector3 momIn; // 3-momentum at entrance to station [GeV/c]
+ TVector3 momOut; // 3-momentum at exit of station [GeV/c]
+ double time = undef::kD; // Time of MC point (after beginning of the event) [ns]
+ int iTmcExt = undef::kI32; // Track ID in the external container
+ // MVD
+ if constexpr (L1DetectorID::kMvd == DetID) {
+ auto* pExtPoint = dynamic_cast<CbmMvdPoint*>(fpMvdPoints->Get(iFile, iEvent, iExtId));
+ if (!pExtPoint) {
+ LOG(warn) << "CbmCaMCModule: MVD MC point with iExtId = " << iExtId << ", iEvent = " << iEvent
+ << ", iFile = " << iFile << " does not exist";
+ return false;
+ }
+ pExtPoint->Position(posIn);
+ pExtPoint->PositionOut(posOut);
+ pExtPoint->Momentum(momIn);
+ pExtPoint->MomentumOut(momOut);
+ time = pExtPoint->GetTime();
+ iTmcExt = pExtPoint->GetTrackID();
+ }
+ // STS
+ else if constexpr (L1DetectorID::kSts == DetID) {
+ auto* pExtPoint = dynamic_cast<CbmStsPoint*>(fpStsPoints->Get(iFile, iEvent, iExtId));
+ if (!pExtPoint) {
+ LOG(warn) << "CbmCaMCModule: STS MC point with iExtId = " << iExtId << ", iEvent = " << iEvent
+ << ", iFile = " << iFile << " does not exist";
+ return false;
+ }
+ pExtPoint->Position(posIn);
+ pExtPoint->PositionOut(posOut);
+ pExtPoint->Momentum(momIn);
+ pExtPoint->MomentumOut(momOut);
+ time = pExtPoint->GetTime();
+ iTmcExt = pExtPoint->GetTrackID();
+ }
+ // MuCh
+ else if constexpr (L1DetectorID::kMuch == DetID) {
+ auto* pExtPoint = dynamic_cast<CbmMuchPoint*>(fpMuchPoints->Get(iFile, iEvent, iExtId));
+ if (!pExtPoint) {
+ LOG(warn) << "CbmCaMCModule: MuCh MC point with iExtId = " << iExtId << ", iEvent = " << iEvent
+ << ", iFile = " << iFile << " does not exist";
+ return false;
+ }
+ pExtPoint->Position(posIn);
+ pExtPoint->PositionOut(posOut);
+ pExtPoint->Momentum(momIn);
+ pExtPoint->Momentum(momOut);
+ time = pExtPoint->GetTime();
+ iTmcExt = pExtPoint->GetTrackID();
+ }
+ // TRD
+ else if constexpr (L1DetectorID::kTrd == DetID) {
+ auto* pExtPoint = dynamic_cast<CbmTrdPoint*>(fpTrdPoints->Get(iFile, iEvent, iExtId));
+ if (!pExtPoint) {
+ LOG(warn) << "CbmCaMCModule: TRD MC point with iExtId = " << iExtId << ", iEvent = " << iEvent
+ << ", iFile = " << iFile << " does not exist";
+ return false;
+ }
+ pExtPoint->Position(posIn);
+ pExtPoint->PositionOut(posOut);
+ pExtPoint->Momentum(momIn);
+ pExtPoint->MomentumOut(momOut);
+ time = pExtPoint->GetTime();
+ iTmcExt = pExtPoint->GetTrackID();
+ }
+ // TOF
+ else if constexpr (L1DetectorID::kTof == DetID) {
+ auto* pExtPoint = dynamic_cast<CbmTofPoint*>(fpTofPoints->Get(iFile, iEvent, iExtId));
+ if (!pExtPoint) {
+ LOG(warn) << "CbmCaMCModule: TOF MC point with iExtId = " << iExtId << ", iEvent = " << iEvent
+ << ", iFile = " << iFile << " does not exist";
+ return false;
+ }
+ pExtPoint->Position(posIn);
+ pExtPoint->Position(posOut);
+ pExtPoint->Momentum(momIn);
+ pExtPoint->Momentum(momOut);
+ time = pExtPoint->GetTime();
+ iTmcExt = pExtPoint->GetTrackID();
+ }
+ if (iTmcExt < 0) {
+ LOG(warn) << "CbmCaMCModule: For MC point with iExtId = " << iExtId << ", iEvent = " << iEvent
+ << ", iFile = " << iFile << " MC track is undefined (has ID = " << iTmcExt << ')';
+ return false;
+ }
+ TVector3 posMid = 0.5 * (posIn + posOut);
+ TVector3 momMid = 0.5 * (momIn + momOut);
+
+ // // ----- Get station index
+ int stationID = undef::kI32; // Index of station in the array of active tracking stations
+ int nStationsGeo = fpParameters->GetNstationsGeometry(DetID);
+ // MVD, STS TODO: Try to extend for MuCh and TRD
+ if constexpr (L1DetectorID::kMvd == DetID || L1DetectorID::kSts == DetID) {
+ // Determine active station global index
+ double bestDist = 1.e+20;
+ for (int iStLocGeo = 0; iStLocGeo < nStationsGeo; ++iStLocGeo) {
+ int iStActive = fpParameters->GetStationIndexActive(iStLocGeo, DetID);
+ if (iStActive < 0) { continue; }
+ const L1Station& station = fpParameters->GetStation(iStActive);
+ double dist = posIn.Z() - station.fZ[0];
+ if (fabs(dist) < fabs(bestDist)) {
+ bestDist = dist;
+ stationID = iStActive;
+ }
+ }
+ }
+ // MuCh, TRD
+ else if constexpr (L1DetectorID::kMuch == DetID || L1DetectorID::kTrd == DetID) {
+ // Offset for MuCh - 2.5 cm, for TRD - 20. cm. TODO: Where did these values came from?
+ double offset = (L1DetectorID::kMuch == DetID) ? 2.5 : 20.;
+ for (int iStLocGeo = 0; iStLocGeo < nStationsGeo; ++iStLocGeo) {
+ int iStActive = fpParameters->GetStationIndexActive(iStLocGeo, DetID);
+ if (iStActive < 0) { continue; }
+ const L1Station& station = fpParameters->GetStation(iStActive);
+ if (posMid.Z() > station.fZ[0] - offset) { stationID = iStActive; }
+ }
+ }
+
+ // Update point time with event time
+ time += fpMCEventList->GetEventTime(iEvent, iFile);
+
+ // ----- Reject MC points falling out of the time slice
+ // STS, MuCh, TRD, TOF
+ if constexpr (DetID != L1DetectorID::kMvd) {
+ double startT = fpTimeSlice->GetStartTime();
+ double endT = fpTimeSlice->GetEndTime();
+
+ if ((startT > 0. && time < startT) || (endT > 0. && time > endT)) {
+ LOG(warn) << "CbmCaMCModule: MC point with iExtId = " << iExtId << ", iEvent = " << iEvent
+ << ", iFile = " << iFile << " and det id " << int(DetID) << " fell out of the TS duration [" << startT
+ << ", " << endT << "] with measured time = " << time << " [ns]";
+ return false;
+ }
+ }
+
+ // ----- Fill MC point
+ point.SetExternalId(iExtId);
+ point.SetEventId(iEvent);
+ point.SetFileId(iFile);
+ point.SetTime(time);
+ point.SetX(posMid.X());
+ point.SetY(posMid.Y());
+ point.SetZ(posMid.Z());
+ point.SetXIn(posIn.X());
+ point.SetYIn(posIn.Y());
+ point.SetZIn(posIn.Z());
+ point.SetXOut(posOut.X());
+ point.SetYOut(posOut.Y());
+ point.SetZOut(posOut.Z());
+ point.SetPx(momMid.X());
+ point.SetPy(momMid.Y());
+ point.SetPz(momMid.Z());
+ point.SetPxIn(momIn.X());
+ point.SetPyIn(momIn.Y());
+ point.SetPzIn(momIn.Z());
+ point.SetPxOut(momOut.X());
+ point.SetPyOut(momOut.Y());
+ point.SetPzOut(momOut.Z());
+
+ int iTmcInt = fpMCData->FindInternalTrackIndex(iTmcExt, iEvent, iFile);
+
+ point.SetId(fpMCData->GetNofPoints()); // select current number of points as a local id of points
+ point.SetTrackId(iTmcInt);
+ point.SetStationId(stationID);
+ point.SetDetectorId(DetID);
+
+ auto* pExtTrk = L1_DYNAMIC_CAST<CbmMCTrack*>(fpMCTracks->Get(iFile, iEvent, iTmcExt));
+ if (!pExtTrk) {
+ LOG(warn) << "CbmCaMCModule: MC track with iTmcExt = " << iTmcExt << ", iEvent = " << iEvent
+ << ", iFile = " << iFile << " MC track is undefined (nullptr)";
+ return false;
+ }
+ point.SetPdgCode(pExtTrk->GetPdgCode());
+ point.SetMotherId(pExtTrk->GetMotherId());
+
+ auto* pPdgDB = TDatabasePDG::Instance()->GetParticle(point.GetPdgCode());
+ point.SetMass(pPdgDB ? pPdgDB->Mass() : 0.); /// TODO: Set from track
+ point.SetCharge(pPdgDB ? pPdgDB->Charge() / 3. : 0.);
+
+ return true;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+// NOTE: template is used, because another template function FillMCPoint is used inside
+template<L1DetectorID DetID>
+void cbm::ca::MCModule::ReadMCPointsForDetector(CbmMCDataArray* pPoints)
+{
+ if constexpr (L1DetectorID::kTof != DetID) {
+ for (const auto& key : fFileEventIDs) {
+ int iFile = key.first;
+ int iEvent = key.second;
+ int nPointsEvent = pPoints->Size(iFile, iEvent);
+ for (int iP = 0; iP < nPointsEvent; ++iP) {
+ ::ca::tools::MCPoint aPoint;
+ if (FillMCPoint<DetID>(iP, iEvent, iFile, aPoint)) {
+ aPoint.SetExternalId(fpMCData->GetPointGlobExtIndex(DetID, iP));
+ int iPInt = fpMCData->GetNofPoints(); // assign an index of point in internal container
+ if (aPoint.GetTrackId() > -1) { fpMCData->GetTrack(aPoint.GetTrackId()).AddPointIndex(iPInt); }
+ fpMCData->AddPoint(aPoint);
+ }
+ } // iP: end
+ } // key: end
+ }
+}
+
+
+#endif // CbmCaMCModule_h
diff --git a/reco/L1/CbmCaTimeSliceReader.cxx b/reco/L1/CbmCaTimeSliceReader.cxx
index 85e0ec20d3b1de721cfedfac9d8e1c97269b983c..bb96dcf479cc181a213a0a8d2d4614853aeca1c7 100644
--- a/reco/L1/CbmCaTimeSliceReader.cxx
+++ b/reco/L1/CbmCaTimeSliceReader.cxx
@@ -19,21 +19,19 @@
#include "FairRootManager.h"
#include "Logger.h"
+#include <algorithm>
#include <numeric>
#include "L1Constants.h"
#include "L1InputData.h"
#include "L1Parameters.h"
+using ca::tools::HitRecord;
using cbm::ca::TimeSliceReader;
using L1Constants::clrs::kCL; // clear log
using L1Constants::clrs::kGNb; // green bold log
using L1Constants::clrs::kRDb; // red bold log
-// ---------------------------------------------------------------------------------------------------------------------
-//
-TimeSliceReader::TimeSliceReader(ECbmTrackingMode mode) : fTrackingMode(mode) {}
-
// ---------------------------------------------------------------------------------------------------------------------
//
void TimeSliceReader::Clear()
@@ -45,11 +43,7 @@ void TimeSliceReader::Clear()
fpBrTimeSlice = nullptr;
fpParameters = nullptr;
- fpBrMvdHits = nullptr;
- fpBrStsHits = nullptr;
- fpBrMuchHits = nullptr;
- fpBrTrdHits = nullptr;
- fpBrTofHits = nullptr;
+ fvpBrHits.fill(nullptr);
fpBrRecoTracks = nullptr;
fpBrStsTracks = nullptr;
@@ -80,17 +74,17 @@ void TimeSliceReader::CheckInit() const
if (!fpBrTimeSlice) { throw std::logic_error("Time slice branch is unavailable"); }
- if (fvbUseDet[L1DetectorID::kMvd] && !fpBrMvdHits) { throw std::logic_error("MVD hits branch is unavailable"); }
- if (fvbUseDet[L1DetectorID::kSts] && !fpBrStsHits) { throw std::logic_error("STS hits branch is unavailable"); }
- if (fvbUseDet[L1DetectorID::kMuch] && !fpBrMuchHits) { throw std::logic_error("MuCh hits branch is unavailable"); }
- if (fvbUseDet[L1DetectorID::kTrd] && !fpBrTrdHits) { throw std::logic_error("TRD hits branch is unavailable"); }
- if (fvbUseDet[L1DetectorID::kTof] && !fpBrTofHits) { throw std::logic_error("TOF hits branch is unavailable"); }
+ for (int iDet = 0; iDet < static_cast<int>(L1DetectorID::kEND); ++iDet) {
+ if (fvbUseDet[iDet] && !fvpBrHits[iDet]) {
+ throw std::logic_error(std::string(kDetName[iDet]) + " hits branch is not found");
+ }
+ }
if (fpvTracks) {
- if (ECbmTrackingMode::kSTS == fTrackingMode) {
+ if (ECbmCaTrackingMode::kSTS == fTrackingMode) {
if (!fpBrRecoTracks) { throw std::logic_error("StsTrack branch is unavailable"); }
}
- else if (ECbmTrackingMode::kMCBM == fTrackingMode) {
+ else if (ECbmCaTrackingMode::kMCBM == fTrackingMode) {
if (!fpBrRecoTracks) { throw std::logic_error("GlobalTrack branch is unavailable"); }
if (fvbUseDet[L1DetectorID::kSts] && !fpBrStsTracks) { throw std::logic_error("StsTrack branch is not found"); }
if (fvbUseDet[L1DetectorID::kMuch] && !fpBrMuchTracks) {
@@ -109,11 +103,11 @@ try {
LOG(info) << "TimeSliceReader: initializing run ... ";
// Init tracking detector interfaces
- fpMvdInterface = CbmMvdTrackingInterface::Instance();
- fpStsInterface = CbmStsTrackingInterface::Instance();
- fpMuchInterface = CbmMuchTrackingInterface::Instance();
- fpTrdInterface = CbmTrdTrackingInterface::Instance();
- fpTofInterface = CbmTofTrackingInterface::Instance();
+ fvpDetInterface[L1DetectorID::kMvd] = CbmMvdTrackingInterface::Instance();
+ fvpDetInterface[L1DetectorID::kSts] = CbmStsTrackingInterface::Instance();
+ fvpDetInterface[L1DetectorID::kMuch] = CbmMuchTrackingInterface::Instance();
+ fvpDetInterface[L1DetectorID::kTrd] = CbmTrdTrackingInterface::Instance();
+ fvpDetInterface[L1DetectorID::kTof] = CbmTofTrackingInterface::Instance();
// ** Init data branches **
@@ -122,22 +116,25 @@ try {
fpBrTimeSlice = dynamic_cast<CbmTimeSlice*>(pFairManager->GetObject("TimeSlice."));
- // Init hit branches
- if (fvbUseDet[L1DetectorID::kMvd]) { fpBrMvdHits = dynamic_cast<TClonesArray*>(pFairManager->GetObject("MvdHit")); }
- if (fvbUseDet[L1DetectorID::kSts]) { fpBrStsHits = dynamic_cast<TClonesArray*>(pFairManager->GetObject("StsHit")); }
- if (fvbUseDet[L1DetectorID::kMuch]) {
- fpBrMuchHits = dynamic_cast<TClonesArray*>(pFairManager->GetObject("MuchPixelHit"));
- }
- if (fvbUseDet[L1DetectorID::kTrd]) { fpBrTrdHits = dynamic_cast<TClonesArray*>(pFairManager->GetObject("TrdHit")); }
- if (fvbUseDet[L1DetectorID::kTof]) { fpBrTofHits = dynamic_cast<TClonesArray*>(pFairManager->GetObject("TofHit")); }
+ // Init branches
+ auto InitHitBranch = [&](L1DetectorID detID, const char* branchName) {
+ if (fvbUseDet[detID]) { fvpBrHits[detID] = dynamic_cast<TClonesArray*>(pFairManager->GetObject(branchName)); }
+ };
+
+ InitHitBranch(L1DetectorID::kMvd, "MvdHit");
+ InitHitBranch(L1DetectorID::kSts, "StsHit");
+ InitHitBranch(L1DetectorID::kMuch, "MuchPixelHit");
+ InitHitBranch(L1DetectorID::kTrd, "TrdHit");
+ InitHitBranch(L1DetectorID::kTof, "TofHit");
// Init track branches
+
if (fpvTracks) {
switch (fTrackingMode) {
- case ECbmTrackingMode::kSTS:
+ case ECbmCaTrackingMode::kSTS:
fpBrRecoTracks = dynamic_cast<TClonesArray*>(pFairManager->GetObject("StsTrack"));
break;
- case ECbmTrackingMode::kMCBM:
+ case ECbmCaTrackingMode::kMCBM:
fpBrRecoTracks = dynamic_cast<TClonesArray*>(pFairManager->GetObject("GlobalTrack"));
if (fvbUseDet[L1DetectorID::kSts]) {
fpBrStsTracks = dynamic_cast<TClonesArray*>(pFairManager->GetObject("StsTrack"));
@@ -188,7 +185,7 @@ void TimeSliceReader::ReadRecoTracks()
int nTracks = fpBrRecoTracks->GetEntriesFast();
fpvTracks->reset(nTracks);
switch (fTrackingMode) {
- case ECbmTrackingMode::kSTS:
+ case ECbmCaTrackingMode::kSTS:
// Fill tracks from StsTrack branch
for (int iT = 0; iT < nTracks; ++iT) {
auto* pInputTrack = static_cast<CbmStsTrack*>(fpBrRecoTracks->At(iT));
@@ -207,18 +204,18 @@ void TimeSliceReader::ReadRecoTracks()
track.Hits.reserve(pInputTrack->GetNofHits());
for (int iH = 0; iH < pInputTrack->GetNofMvdHits(); ++iH) {
int iHext = pInputTrack->GetMvdHitIndex(iH);
- int iHint = vHitMvdIds[iHext];
+ int iHint = fvvHitExtToIntIndexMap[L1DetectorID::kMvd][iHext];
track.Hits.push_back(iHint);
} // iH
for (int iH = 0; iH < pInputTrack->GetNofStsHits(); ++iH) {
int iHext = pInputTrack->GetStsHitIndex(iH);
- int iHint = vHitStsIds[iHext];
+ int iHint = fvvHitExtToIntIndexMap[L1DetectorID::kSts][iHext];
track.Hits.push_back(iHint);
} // iH
} // iT
break;
- case ECbmTrackingMode::kMCBM:
+ case ECbmCaTrackingMode::kMCBM:
//LOG(fatal) << "Sorry, mCBM mode has not been implemented yet. Stay tuned :)";
// Fill tracks from GlobalTrack branch
for (int iT = 0; iT < nTracks; ++iT) {
@@ -238,13 +235,13 @@ void TimeSliceReader::ReadRecoTracks()
if (fvbUseDet[L1DetectorID::kMvd]) {
for (int iH = 0; iH < pStsTrack->GetNofMvdHits(); ++iH) {
int iHext = pStsTrack->GetMvdHitIndex(iH);
- int iHint = vHitMvdIds[iHext];
+ int iHint = fvvHitExtToIntIndexMap[L1DetectorID::kMvd][iHext];
track.Hits.push_back(iHint);
}
}
for (int iH = 0; iH < pStsTrack->GetNofStsHits(); ++iH) {
int iHext = pStsTrack->GetStsHitIndex(iH);
- int iHint = vHitStsIds[iHext];
+ int iHint = fvvHitExtToIntIndexMap[L1DetectorID::kSts][iHext];
track.Hits.push_back(iHint);
}
}
@@ -257,7 +254,7 @@ void TimeSliceReader::ReadRecoTracks()
auto* pMuchTrack = static_cast<CbmMuchTrack*>(fpBrMuchTracks->At(iMuchTrkId));
for (int iH = 0; iH < pMuchTrack->GetNofHits(); ++iH) {
int iHext = pMuchTrack->GetHitIndex(iH);
- int iHint = vHitMuchIds[iHext];
+ int iHint = fvvHitExtToIntIndexMap[L1DetectorID::kMuch][iHext];
track.Hits.push_back(iHint);
}
}
@@ -270,7 +267,7 @@ void TimeSliceReader::ReadRecoTracks()
const auto* pTrdTrack = static_cast<const CbmTrdTrack*>(fpBrTrdTracks->At(iTrdTrkId));
for (int iH = 0; iH < pTrdTrack->GetNofHits(); ++iH) {
int iHext = pTrdTrack->GetHitIndex(iH);
- int iHint = vHitTrdIds[iHext];
+ int iHint = fvvHitExtToIntIndexMap[L1DetectorID::kTrd][iHext];
track.Hits.push_back(iHint);
} // iH
}
@@ -283,7 +280,7 @@ void TimeSliceReader::ReadRecoTracks()
const auto* pTofTrack = static_cast<const CbmTofTrack*>(fpBrTofTracks->At(iTofTrkId));
for (int iH = 0; iH < pTofTrack->GetNofHits(); ++iH) {
int iHext = pTofTrack->GetHitIndex(iH);
- int iHint = vHitTofIds[iHext];
+ int iHint = fvvHitExtToIntIndexMap[L1DetectorID::kTof][iHext];
track.Hits.push_back(iHint);
} // iH
} // if iTofTrkId > -1
@@ -303,42 +300,28 @@ void TimeSliceReader::RegisterIODataManager(std::shared_ptr<L1IODataManager>& pI
// ---------------------------------------------------------------------------------------------------------------------
//
-template<>
-void TimeSliceReader::SortQaHits<true>()
+void TimeSliceReader::SortQaHits()
{
int nStationsActive = fpParameters->GetNstationsActive();
L1Vector<CbmL1HitDebugInfo> vNewHits(fpvQaHits->size());
- std::vector<int> vHitFirstIndexes(nStationsActive + 1, 0);
+ std::vector<int> vHitFstIndexes(nStationsActive + 1, 0);
std::vector<int> vNofHitsStored(nStationsActive, 0);
// Count number of hits in each station (NOTE: we could use here boarders from the IO data manager, but we would keep
// these two procedures independent)
- for (const auto& hit : (*fpvQaHits)) {
- ++vHitFirstIndexes[hit.GetStationId() + 1];
- }
-
+ std::for_each(fpvQaHits->begin(), fpvQaHits->end(), [&](const auto& h) { ++vHitFstIndexes[h.GetStationId() + 1]; });
for (int iSt = 0; iSt < nStationsActive; ++iSt) {
- vHitFirstIndexes[iSt + 1] += vHitFirstIndexes[iSt];
+ vHitFstIndexes[iSt + 1] += vHitFstIndexes[iSt];
}
for (const auto& hit : (*fpvQaHits)) {
- int iSt = hit.GetStationId();
- vNewHits[vHitFirstIndexes[iSt] + (vNofHitsStored[iSt]++)] = hit;
+ int iSt = hit.GetStationId();
+ vNewHits[vHitFstIndexes[iSt] + (vNofHitsStored[iSt]++)] = hit;
}
std::swap(vNewHits, (*fpvQaHits));
}
-// ---------------------------------------------------------------------------------------------------------------------
-//
-template<>
-void TimeSliceReader::SortQaHits<false>()
-{
- std::sort(fpvQaHits->begin(), fpvQaHits->end(), [](const CbmL1HitDebugInfo& a, const CbmL1HitDebugInfo& b) {
- return (a.iStation < b.iStation) || (a.iStation == b.iStation) && (a.y < b.y);
- });
-}
-
// ---------------------------------------------------------------------------------------------------------------------
//
void TimeSliceReader::ReadHits()
@@ -347,11 +330,10 @@ void TimeSliceReader::ReadHits()
fNofHitKeys = 0;
fFirstHitKey = 0;
- if (fvbUseDet[L1DetectorID::kMvd]) { fvNofHitsTotal[L1DetectorID::kMvd] = fpBrMvdHits->GetEntriesFast(); }
- if (fvbUseDet[L1DetectorID::kSts]) { fvNofHitsTotal[L1DetectorID::kSts] = fpBrStsHits->GetEntriesFast(); }
- if (fvbUseDet[L1DetectorID::kMuch]) { fvNofHitsTotal[L1DetectorID::kMuch] = fpBrMuchHits->GetEntriesFast(); }
- if (fvbUseDet[L1DetectorID::kTrd]) { fvNofHitsTotal[L1DetectorID::kTrd] = fpBrTrdHits->GetEntriesFast(); }
- if (fvbUseDet[L1DetectorID::kTof]) { fvNofHitsTotal[L1DetectorID::kTof] = fpBrTofHits->GetEntriesFast(); }
+ // TODO: Address case with CbmEvent != nullptr
+ for (int iDet = 0; iDet < static_cast<int>(L1DetectorID::kEND); ++iDet) {
+ if (fvbUseDet[iDet]) { fvNofHitsTotal[iDet] = fvpBrHits[iDet]->GetEntriesFast(); }
+ }
int nHitsTot = std::accumulate(fvNofHitsTotal.begin(), fvNofHitsTotal.end(), 0);
@@ -364,29 +346,16 @@ void TimeSliceReader::ReadHits()
fpvQaHits->clear();
fpvQaHits->reserve(nHitsTot);
}
- if (fpIODataManager) {
- fpIODataManager->ResetInputData();
- fpIODataManager->ReserveNhits(nHitsTot);
- }
+ if (fpIODataManager) { fpIODataManager->ResetInputData(nHitsTot); }
std::fill(fvHitFirstIndexDet.begin(), fvHitFirstIndexDet.end(), 0);
// Read hits for different detectors
- if (fvbUseDet[L1DetectorID::kMvd]) {
- fvNofHitsUsed[L1DetectorID::kMvd] = ReadHitsForDetector<L1DetectorID::kMvd>(fpBrMvdHits);
- }
- if (fvbUseDet[L1DetectorID::kSts]) {
- fvNofHitsUsed[L1DetectorID::kSts] = ReadHitsForDetector<L1DetectorID::kSts>(fpBrStsHits);
- }
- if (fvbUseDet[L1DetectorID::kMuch]) {
- fvNofHitsUsed[L1DetectorID::kMuch] = ReadHitsForDetector<L1DetectorID::kMuch>(fpBrMuchHits);
- }
- if (fvbUseDet[L1DetectorID::kTrd]) {
- fvNofHitsUsed[L1DetectorID::kTrd] = ReadHitsForDetector<L1DetectorID::kTrd>(fpBrTrdHits);
- }
- if (fvbUseDet[L1DetectorID::kTof]) {
- fvNofHitsUsed[L1DetectorID::kTof] = ReadHitsForDetector<L1DetectorID::kTof>(fpBrTofHits);
- }
+ fvNofHitsUsed[L1DetectorID::kMvd] = ReadHitsForDetector<L1DetectorID::kMvd>();
+ fvNofHitsUsed[L1DetectorID::kSts] = ReadHitsForDetector<L1DetectorID::kSts>();
+ fvNofHitsUsed[L1DetectorID::kMuch] = ReadHitsForDetector<L1DetectorID::kMuch>();
+ fvNofHitsUsed[L1DetectorID::kTrd] = ReadHitsForDetector<L1DetectorID::kTrd>();
+ fvNofHitsUsed[L1DetectorID::kTof] = ReadHitsForDetector<L1DetectorID::kTof>();
// Save first hit index for different detector subsystems
for (uint32_t iDet = 0; iDet < fvNofHitsUsed.size(); ++iDet) {
@@ -399,26 +368,18 @@ void TimeSliceReader::ReadHits()
if (fpIODataManager) { fpIODataManager->SetNhitKeys(fNofHitKeys); }
// Sort debug hits
- if (fpvQaHits) { this->SortQaHits<true>(); } // false - old sort, true - new backet sort
+ if (fpvQaHits) { this->SortQaHits(); }
// Update maps of ext->int hit indexes
// NOTE: fvpHitIds must be initialized, if we want to read tracks from the file
if (fpvHitIds) {
- auto UpdateHitIndexMap = [&](L1Vector<int> & vIds, L1DetectorID detID) constexpr
- {
- if (fvbUseDet[detID]) {
- vIds.reset(fvNofHitsTotal[detID]);
- for (int iH = fvHitFirstIndexDet[int(detID)]; iH < fvHitFirstIndexDet[int(detID) + 1]; ++iH) {
- vIds[(*fpvQaHits)[iH].ExtIndex] = iH;
- }
- }
- };
+ auto ResetIndexMap = [&, idx = 0](auto& v) mutable { v.reset(fvNofHitsTotal[idx++]); };
+ std::for_each(fvvHitExtToIntIndexMap.begin(), fvvHitExtToIntIndexMap.end(), ResetIndexMap);
- UpdateHitIndexMap(vHitMvdIds, L1DetectorID::kMvd);
- UpdateHitIndexMap(vHitStsIds, L1DetectorID::kSts);
- UpdateHitIndexMap(vHitMuchIds, L1DetectorID::kMuch);
- UpdateHitIndexMap(vHitTrdIds, L1DetectorID::kTrd);
- UpdateHitIndexMap(vHitTofIds, L1DetectorID::kTof);
+ for (int iH = 0; iH < fNofHits; ++iH) {
+ const auto& hit = (*fpvQaHits)[iH];
+ fvvHitExtToIntIndexMap[hit.Det][hit.ExtIndex] = iH;
+ }
}
}
@@ -448,17 +409,17 @@ void TimeSliceReader::StoreHitRecord(const HitRecord& hitRecord)
// Save debug information
if (fpvQaHits) {
- CbmL1HitDebugInfo hitDbg;
- hitDbg.Det = hitRecord.fDet;
- hitDbg.ExtIndex = hitRecord.fExtId;
- hitDbg.IntIndex = static_cast<int>(fpvQaHits->size());
- hitDbg.iStation = hitRecord.fStaId;
- hitDbg.x = hitRecord.fX;
- hitDbg.y = hitRecord.fY;
- hitDbg.dx = hitRecord.fDx;
- hitDbg.dy = hitRecord.fDy;
- hitDbg.dxy = hitRecord.fDxy;
- hitDbg.time = hitRecord.fT;
- fpvQaHits->push_back(hitDbg);
+ CbmL1HitDebugInfo aHitQa;
+ aHitQa.Det = hitRecord.fDet;
+ aHitQa.ExtIndex = hitRecord.fExtId;
+ aHitQa.IntIndex = static_cast<int>(fpvQaHits->size());
+ aHitQa.iStation = hitRecord.fStaId;
+ aHitQa.x = hitRecord.fX;
+ aHitQa.y = hitRecord.fY;
+ aHitQa.dx = hitRecord.fDx;
+ aHitQa.dy = hitRecord.fDy;
+ aHitQa.dxy = hitRecord.fDxy;
+ aHitQa.time = hitRecord.fT;
+ fpvQaHits->push_back(aHitQa);
}
}
diff --git a/reco/L1/CbmCaTimeSliceReader.h b/reco/L1/CbmCaTimeSliceReader.h
index 1b9e60dd05fc69b666c191301ecd081fe96d72c8..7cab89abfe338c071b108291cf683ba39cffeaba 100644
--- a/reco/L1/CbmCaTimeSliceReader.h
+++ b/reco/L1/CbmCaTimeSliceReader.h
@@ -10,6 +10,7 @@
#ifndef CbmCaTimeSliceReader_h
#define CbmCaTimeSliceReader_h 1
+#include "CbmCaMCModule.h"
#include "CbmL1.h" // TMP: for CbmL1HitStore
#include "CbmL1DetectorID.h"
#include "CbmL1Hit.h"
@@ -28,6 +29,7 @@
#include "TClonesArray.h"
+#include "CaToolsHitRecord.h"
#include "L1Constants.h"
#include "L1Vector.h"
@@ -38,39 +40,19 @@ class L1Parameters;
namespace cbm::ca
{
+ using ::ca::tools::HitRecord;
+
+
/// @brief A reader of time slice for CA tracker
///
/// The class reads reconstructed hits and reconstructed tracks (optionally) and fills the CA tracking internal
/// data structures.
///
class TimeSliceReader {
- /// @brief A structure to store hits information from different detectors in a uniform manner
- ///
- struct HitRecord {
- double fX = 0.; ///< x component of hit position [cm]
- double fY = 0.; ///< y component of hit position [cm]
- double fDx = 0.; ///< error of x component of hit position [cm]
- double fDy = 0.; ///< error of y component of hit position [cm]
- double fDxy = 0.; ///< correlation between x and y components [cm]
- double fU = 0.; ///< hit position in direction of front strips [cm]
- double fV = 0.; ///< hit position in direction of back strips [cm]
- double fDu = 0.; ///< hit position error in direction of front strips [cm]
- double fDv = 0.; ///< hit position error in direction of back strips [cm]
- double fZ = 0.; ///< z component of hit position [cm]
- double fT = 0.; ///< time of hit [ns]
- double fDt = 0.; ///< time error of hit [ns]
- int64_t fDataStream = -1; ///< Global index of detector module
- int fExtId = -1; ///< external index of hit
- int fStaId = -1; ///< index of active tracking station
- int fStripF = -1; ///< index of front strip
- int fStripB = -1; ///< index of back strip
- int fDet = -1; ///< detector ID
- };
-
public:
/// @brief Constructor from parameters
/// @param mode Tracking mode
- TimeSliceReader(ECbmTrackingMode mode);
+ TimeSliceReader() = default;
/// @brief Destructor
~TimeSliceReader() = default;
@@ -90,10 +72,6 @@ namespace cbm::ca
/// @brief Clears class content
void Clear();
- /// @brief Check class initialization
- /// @note The function throws std::logic_error, if initialization is incomplete
- void CheckInit() const;
-
/// @brief Gets reference to container of first hit indexes in a detector subsystem
/// @return Ref. to the container
const auto& GetHitFirstIndexDet() const { return fvHitFirstIndexDet; }
@@ -106,6 +84,9 @@ namespace cbm::ca
return fvHitFirstIndexDet[int(iDet) + 1] - fvHitFirstIndexDet[int(iDet)];
}
+ /// @brief Gets CBM tracking mode
+ ECbmCaTrackingMode GetTrackingMode() const { return fTrackingMode; }
+
/// @brief Run initializer function
/// @return Success flag
///
@@ -117,12 +98,6 @@ namespace cbm::ca
/// Reads hits and tracks (optionally) from time slice
void InitTimeSlice();
- /// @brief Reads hits
- void ReadHits();
-
- /// @brief Reads reconstructed tracks
- void ReadRecoTracks();
-
/// @brief Registers hit debug info container
/// @param vQaHits Reference to Qa hit container
/// @note If no container is registered, all related routines are omitted
@@ -153,17 +128,30 @@ namespace cbm::ca
/// @note Should be called before this->Init()
void SetDetector(L1DetectorID detID, bool flag = true) { fvbUseDet[detID] = flag; }
+ /// @brief Sets the tracking mode
+ /// @param mode Tracking mode (from ECbmTrackingMode)
+ void SetTrackingMode(ECbmCaTrackingMode mode) { fTrackingMode = mode; }
+
+
private:
+ /// @brief Check class initialization
+ /// @note The function throws std::logic_error, if initialization is incomplete
+ void CheckInit() const;
+
+ /// @brief Reads hits
+ void ReadHits();
+
+ /// @brief Reads reconstructed tracks
+ void ReadRecoTracks();
+
/// @brief Reads hits for a given detector subsystem
/// @tparam Detector ID
- /// @param pBrHits Pointer to input hit branch
/// @return Number of stored hits
/// @note The function modifies fNofHitKey and fFirstHitKey counters
template<L1DetectorID DetID>
- int ReadHitsForDetector(const TClonesArray* pBrHits);
+ int ReadHitsForDetector();
/// @brief Sorts QA hit objects by stations
- template<bool IsNewSortingApproach>
void SortQaHits();
/// @brief Saves hit to data structures
@@ -172,27 +160,14 @@ namespace cbm::ca
/// Stores recorded hit information into registered hit containers
void StoreHitRecord(const HitRecord& hitRecord);
-
bool fbReadTracks = true; ///< flag to read reconstructed tracks from reco.root
-
- // Pointers to the tracking detector interfaces
- CbmTrackingDetectorInterfaceBase* fpMvdInterface = nullptr;
- CbmTrackingDetectorInterfaceBase* fpStsInterface = nullptr;
- CbmTrackingDetectorInterfaceBase* fpMuchInterface = nullptr;
- CbmTrackingDetectorInterfaceBase* fpTrdInterface = nullptr;
- CbmTrackingDetectorInterfaceBase* fpTofInterface = nullptr;
+ /// @brief Pointers to the tracking detector interfaces for each subsystem
+ CbmCaDetIdArr_t<const CbmTrackingDetectorInterfaceBase*> fvpDetInterface = {{nullptr}};
// Input data branches
- CbmTimeSlice* fpBrTimeSlice = nullptr; ///< Pointer to the TS object
-
- TClonesArray* fpBrMvdHits = nullptr; ///< Input branch for MVD hits ("MvdHit")
- TClonesArray* fpBrStsHits = nullptr; ///< Input branch for STS hits ("StsHit")
- TClonesArray* fpBrMuchHits = nullptr; ///< Input branch for MuCh hits ("MuchPixelHit")
- TClonesArray* fpBrTrdHits = nullptr; ///< Input branch for TRD hits ("TrdHit")
- TClonesArray* fpBrTofHits = nullptr; ///< Input branch for TOF hits ("TofHit")
-
- CbmCaDetIdArr_t<TClonesArray*> fvpBrHits = {nullptr}; ///< Input branch for hits
+ CbmTimeSlice* fpBrTimeSlice = nullptr; ///< Pointer to the TS object
+ CbmCaDetIdArr_t<TClonesArray*> fvpBrHits = {{nullptr}}; ///< Input branch for hits
// Branches for reconstructed tracks. The input at the moment (as for 27.02.2023) depends on the selected
// tracking mode. For simulations in CBM, the CA tracking is used only in STS + MVD detectors. In this case
@@ -211,27 +186,21 @@ namespace cbm::ca
std::shared_ptr<L1IODataManager> fpIODataManager = nullptr; ///< Pointer to input data manager
std::shared_ptr<L1Parameters> fpParameters = nullptr; ///< Pointer to tracking parameters object
-
// Maps of hit indexes: ext -> int
- L1Vector<int> vHitMvdIds {"CbmCaTimeSliceReader::vHitMvdIds"};
- L1Vector<int> vHitStsIds {"CbmCaTimeSliceReader::vHitStsIds"};
- L1Vector<int> vHitMuchIds {"CbmCaTimeSliceReader::vHitMuchIds"};
- L1Vector<int> vHitTrdIds {"CbmCaTimeSliceReader::vHitTrdIds"};
- L1Vector<int> vHitTofIds {"CbmCaTimeSliceReader::vHitTofIds"};
+ CbmCaDetIdArr_t<L1Vector<int>> fvvHitExtToIntIndexMap; ///< Hit index map ext -> int
- CbmCaDetIdArr_t<int> fvNofHitsTotal = {0}; ///< Total hit number in detector
- CbmCaDetIdArr_t<int> fvNofHitsUsed = {0}; ///< Number of used hits in detector
- CbmCaDetIdArr_t<bool> fvbUseDet = {false}; ///< Flag: is detector subsystem used
+ CbmCaDetIdArr_t<int> fvNofHitsTotal = {{0}}; ///< Total hit number in detector
+ CbmCaDetIdArr_t<int> fvNofHitsUsed = {{0}}; ///< Number of used hits in detector
+ CbmCaDetIdArr_t<bool> fvbUseDet = {{false}}; ///< Flag: is detector subsystem used
- // Additional
- ECbmTrackingMode fTrackingMode; ///< Tracking mode
+ ECbmCaTrackingMode fTrackingMode = ECbmCaTrackingMode::kSTS; ///< Tracking mode
// Variables for storing cache
int fNofHits = 0; ///< Stored number of hits
int fNofHitKeys = 0; ///< Recorded number of hit keys
int fFirstHitKey = 0; ///< First index of hit key for the detector subsystem
- std::array<int, L1Constants::size::kMaxNdetectors + 1> fvHitFirstIndexDet = {0}; ///< First hit index in detector
+ std::array<int, L1Constants::size::kMaxNdetectors + 1> fvHitFirstIndexDet = {{0}}; ///< First hit index in detector
};
} // namespace cbm::ca
@@ -243,78 +212,61 @@ namespace cbm::ca
// ---------------------------------------------------------------------------------------------------------------------
//
template<L1DetectorID DetID>
-int cbm::ca::TimeSliceReader::ReadHitsForDetector(const TClonesArray* pBrHits)
+int cbm::ca::TimeSliceReader::ReadHitsForDetector()
{
- int nHitsTot = pBrHits->GetEntriesFast(); // total number of hits stored in a branch
- int nHitsStored = 0; // number of hits stored
+ if (!fvbUseDet[DetID]) { return 0; } // Detector is entirelly not used
+
+ const auto* pDetInterface = fvpDetInterface[DetID];
+ int nHitsTot = fvpBrHits[DetID]->GetEntriesFast(); // total number of hits provided by hit finder
+ int nHitsStored = 0; // number of hits used in tracking
fFirstHitKey = fNofHitKeys;
for (int iHext = 0; iHext < nHitsTot; ++iHext) {
- HitRecord hitRecord; // Record of hits information
+ HitRecord hitRecord;
CbmPixelHit* pPixelHit = nullptr; // Pointer to hit object
- float phiF = 0.; // Stereo angle of front strips
- float phiB = 0.; // Stereo angle of back strips
int iStGeom = -1; // Geometry station number
// Fill out detector specific data
if constexpr (L1DetectorID::kMvd == DetID) {
- CbmMvdHit* pMvdHit = static_cast<CbmMvdHit*>(pBrHits->At(iHext));
- iStGeom = fpMvdInterface->GetTrackingStationIndex(pMvdHit->GetStationNr());
- phiF = fpMvdInterface->GetStripsStereoAngleFront(iStGeom);
- phiB = fpMvdInterface->GetStripsStereoAngleBack(iStGeom);
+ CbmMvdHit* pMvdHit = static_cast<CbmMvdHit*>(fvpBrHits[DetID]->At(iHext));
pPixelHit = static_cast<CbmPixelHit*>(pMvdHit);
- hitRecord.fDu = pMvdHit->GetDx();
- hitRecord.fDv = pMvdHit->GetDy();
+ iStGeom = pDetInterface->GetTrackingStationIndex(pMvdHit->GetStationNr());
}
else if constexpr (L1DetectorID::kSts == DetID) {
- CbmStsHit* pStsHit = static_cast<CbmStsHit*>(pBrHits->At(iHext));
- iStGeom = fpStsInterface->GetTrackingStationIndex(pStsHit->GetAddress());
- phiF = fpStsInterface->GetStripsStereoAngleFront(iStGeom);
- phiB = fpStsInterface->GetStripsStereoAngleBack(iStGeom);
+ CbmStsHit* pStsHit = static_cast<CbmStsHit*>(fvpBrHits[DetID]->At(iHext));
pPixelHit = static_cast<CbmPixelHit*>(pStsHit);
+ iStGeom = pDetInterface->GetTrackingStationIndex(pStsHit->GetAddress());
hitRecord.fStripF = fFirstHitKey + pStsHit->GetFrontClusterId();
hitRecord.fStripB = fFirstHitKey + pStsHit->GetBackClusterId();
hitRecord.fDu = pStsHit->GetDu();
hitRecord.fDv = pStsHit->GetDv();
}
else if constexpr (L1DetectorID::kMuch == DetID) {
- CbmMuchPixelHit* pMuchHit = static_cast<CbmMuchPixelHit*>(pBrHits->At(iHext));
- iStGeom = fpMuchInterface->GetTrackingStationIndex(pMuchHit->GetAddress());
- phiF = fpMuchInterface->GetStripsStereoAngleFront(iStGeom);
- phiB = fpMuchInterface->GetStripsStereoAngleBack(iStGeom);
+ CbmMuchPixelHit* pMuchHit = static_cast<CbmMuchPixelHit*>(fvpBrHits[DetID]->At(iHext));
pPixelHit = static_cast<CbmPixelHit*>(pMuchHit);
- hitRecord.fDu = pMuchHit->GetDx();
- hitRecord.fDv = pMuchHit->GetDy();
+ iStGeom = pDetInterface->GetTrackingStationIndex(pMuchHit->GetAddress());
}
else if constexpr (L1DetectorID::kTrd == DetID) {
- CbmTrdHit* pTrdHit = static_cast<CbmTrdHit*>(pBrHits->At(iHext));
- iStGeom = fpTrdInterface->GetTrackingStationIndex(pTrdHit->GetAddress());
- phiF = fpTrdInterface->GetStripsStereoAngleFront(iStGeom);
- phiB = fpTrdInterface->GetStripsStereoAngleBack(iStGeom);
+ CbmTrdHit* pTrdHit = static_cast<CbmTrdHit*>(fvpBrHits[DetID]->At(iHext));
pPixelHit = static_cast<CbmPixelHit*>(pTrdHit);
- hitRecord.fDu = pTrdHit->GetDx();
- hitRecord.fDv = pTrdHit->GetDy();
+ iStGeom = pDetInterface->GetTrackingStationIndex(pTrdHit->GetAddress());
}
else if constexpr (L1DetectorID::kTof == DetID) {
- CbmTofHit* pTofHit = static_cast<CbmTofHit*>(pBrHits->At(iHext));
+ CbmTofHit* pTofHit = static_cast<CbmTofHit*>(fvpBrHits[DetID]->At(iHext));
+ pPixelHit = static_cast<CbmPixelHit*>(pTofHit);
// NOTE: In TOF we can take station index only from hit, because the function needs information on x and z
// of the hit in case of "beam_mcbm_2021_07_surveyed" (missingHits flag = true).
// TODO: Investigate this case or apply the hack to the TOF level
- iStGeom = fpTofInterface->GetTrackingStationIndex(pTofHit);
- phiF = fpTofInterface->GetStripsStereoAngleFront(iStGeom);
- phiB = fpTofInterface->GetStripsStereoAngleBack(iStGeom);
- pPixelHit = static_cast<CbmPixelHit*>(pTofHit);
- hitRecord.fDu = pTofHit->GetDx();
- hitRecord.fDv = pTofHit->GetDy();
+ iStGeom = pDetInterface->GetTrackingStationIndex(pTofHit);
// *** Additional cuts for TOF ***
// Skip T0 hits
if (5 == CbmTofAddress::GetSmType(pTofHit->GetAddress())) { continue; }
// FIXME: Figure it out, if this cut is still needed (introduced a year ago for mCBM)
- if (ECbmTrackingMode::kMCBM == fTrackingMode && pTofHit->GetZ() > 400) { continue; }
+ if (ECbmCaTrackingMode::kMCBM == fTrackingMode && pTofHit->GetZ() > 400) { continue; }
}
int iStActive = fpParameters->GetStationIndexActive(iStGeom, DetID);
@@ -332,6 +284,8 @@ int cbm::ca::TimeSliceReader::ReadHitsForDetector(const TClonesArray* pBrHits)
hitRecord.fDt = pPixelHit->GetTimeError();
hitRecord.fDet = static_cast<int>(DetID);
hitRecord.fDataStream = (static_cast<int64_t>(hitRecord.fDet) << 60) | pPixelHit->GetAddress();
+ float phiF = pDetInterface->GetStripsStereoAngleFront(iStGeom);
+ float phiB = pDetInterface->GetStripsStereoAngleBack(iStGeom);
hitRecord.fU = hitRecord.fX * cos(phiF) + hitRecord.fY * sin(phiF);
hitRecord.fV = hitRecord.fX * cos(phiB) + hitRecord.fY * sin(phiB);
hitRecord.fExtId = iHext;
@@ -342,6 +296,8 @@ int cbm::ca::TimeSliceReader::ReadHitsForDetector(const TClonesArray* pBrHits)
if (fNofHitKeys <= hitRecord.fStripB) { fNofHitKeys += hitRecord.fStripB; }
}
else {
+ hitRecord.fDu = hitRecord.fDx;
+ hitRecord.fDv = hitRecord.fDy;
hitRecord.fStripF = fFirstHitKey + iHext;
hitRecord.fStripB = hitRecord.fStripF;
if (fNofHitKeys <= hitRecord.fStripF) { fNofHitKeys += hitRecord.fStripF; }
@@ -351,8 +307,8 @@ int cbm::ca::TimeSliceReader::ReadHitsForDetector(const TClonesArray* pBrHits)
this->StoreHitRecord(hitRecord);
++nHitsStored;
} // iH
+
return nHitsStored;
}
-
#endif // CbmCaTimeSliceReader_h
diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx
index d2469ffd5e2e5a7c126b9e92f16ad7d122292fda..0cc4dba49d07b25b681566e4cc0fd4d7f36ab490 100644
--- a/reco/L1/CbmL1.cxx
+++ b/reco/L1/CbmL1.cxx
@@ -1,6 +1,6 @@
/* Copyright (C) 2006-2021 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
SPDX-License-Identifier: GPL-3.0-only
- Authors: Ivan Kisel, Sergey Gorbunov [committer], Valentina Akishina, Igor Kulakov, Maksym Zyzak */
+ Authors: Ivan Kisel, Sergey Gorbunov [committer], Valentina Akishina, Igor Kulakov, Maksym Zyzak, Sergei Zharko */
/*
*====================================================================
@@ -1138,8 +1138,9 @@ void CbmL1::IdealTrackFinder()
for (unsigned int iH = 0; iH < MC.Hits.size(); iH++) {
const int hitI = MC.Hits[iH];
const CbmL1HitDebugInfo& hit = fvHitDebugInfo[hitI];
-
- const int iStation = fvMCPoints[hit.mcPointIds[0]].iStation;
+ const int iP = hit.GetMCPointIndex();
+ assert(iP >= 0);
+ const int iStation = fvMCPoints[iP].iStation;
if (iStation >= 0) hitIndices[iStation] = hitI;
}
diff --git a/reco/L1/CbmL1.h b/reco/L1/CbmL1.h
index 85608d727df3f472967998948eb93713a2009a6b..8d663caccc1343ac81a18ef01f1f5fec3472bd26 100644
--- a/reco/L1/CbmL1.h
+++ b/reco/L1/CbmL1.h
@@ -381,6 +381,9 @@ private:
// ** Reconstruction Performance **
// ********************************
+ /// @brief Sets random seed to CA internal random generator
+ /// @param seed Random seed
+ void SetRandomSeed(unsigned int seed) { fInitManager.SetRandomSeed(seed); }
/// Matches reconstructed and MC tracks
/// \note Should be called before Performances
diff --git a/reco/L1/CbmL1DetectorID.h b/reco/L1/CbmL1DetectorID.h
index 1ec1ae7eb2da97b2dd9bb7f2fd6c1665d271fffd..d4435388cc21589f74802c741838c0c5edf2b829 100644
--- a/reco/L1/CbmL1DetectorID.h
+++ b/reco/L1/CbmL1DetectorID.h
@@ -10,12 +10,15 @@
#ifndef CbmL1DetectorID_h
#define CbmL1DetectorID_h 1
+#include <string>
+
#include "L1EArray.h"
/// Enumeration for the detector subsystems used in L1 tracking
/// It is important for the subsystems to be specified in the actual order. The order is used
/// for the L1Station array filling.
/// Note: L1DetectorID has a forward declaration in L1InitManager.h and L1BaseStationInfo.h
+/// TODO: Probably, we need to move everything into a single CbmCaConst.h header
enum class L1DetectorID
{
kMvd,
@@ -31,9 +34,15 @@ enum class L1DetectorID
template<typename T>
using CbmCaDetIdArr_t = L1EArray<L1DetectorID, T>;
+namespace cbm::ca
+{
+ /// @brief Name of detector subsystems
+ constexpr CbmCaDetIdArr_t<const char*> kDetName = {{"MVD", "STS", "MuCh", "TRD", "TOF"}};
+
+} // namespace cbm::ca
/// @brief Enumeration for different tracking running modes
-enum class ECbmTrackingMode
+enum class ECbmCaTrackingMode
{
kSTS, ///< Local tracking in CBM (STS + MVD), results stored to the StsTrack branch
kMCBM ///< Global tracking in mCBM (STS, MuCh, TRD, TOF), results stored to GlobalTrack branch
diff --git a/reco/L1/CbmL1Hit.h b/reco/L1/CbmL1Hit.h
index 3c5c66caf91ebd9da4569273dbe1f7e2a9ec3d99..92b76f8770f3458d57cb2c7a111c84f007639396 100644
--- a/reco/L1/CbmL1Hit.h
+++ b/reco/L1/CbmL1Hit.h
@@ -49,7 +49,7 @@ public:
///
class CbmL1HitDebugInfo { // TODO: SZh 21.09.2022: Replace instances of this class with L1Hit
public:
- /// Gets detector type
+ /// @brief Gets detector type
/// 0 - MVD
/// 1 - STS
/// 2 - MuCh
@@ -58,29 +58,31 @@ public:
// TODO: SZh 02.03.2023: Replace det. ID with L1DetectorID
int GetDetectorType() const { return Det; }
- /// Gets distance from z-axis [cm]
- float GetR() const { return std::sqrt(x * x + y * y); }
+ /// @brief Gets distance from z-axis [cm]
+ double GetR() const { return std::sqrt(x * x + y * y); }
- /// Gets index of the hit in the external container
+ /// @brief Gets index of the hit in the external container
int GetExternalHitId() const { return ExtIndex; }
- /// Gets MC point index container
- const auto& GetMCPointIndexes() const { return mcPointIds; }
+ /// @brief Gets index of matched MC point
+ int GetMCPointIndex() const { return fPointID; }
- /// Gets index of matched MC point
- int GetMCPointIndex() const { return mcPointIds.size() ? mcPointIds[0] : -1; }
-
- /// Gets global index of active tracking station
+ /// @brief Gets global index of active tracking station
int GetStationId() const { return iStation; }
- /// Gets time measurement [ns]
- float GetT() const { return time; }
+ /// @brief Gets time measurement [ns]
+ double GetT() const { return time; }
+
+ /// @brief Gets x component of position [cm]
+ double GetX() const { return x; }
+
+ /// @brief Gets y component of position [cm]
+ double GetY() const { return y; }
- /// Gets x component of position [cm]
- float GetX() const { return x; }
+ /// @brief Sets index of matched MC point
+ /// @param pointID Internal index of MC point
+ void SetMCPointIndex(int pointID) { fPointID = pointID; }
- /// Gets y component of position [cm]
- float GetY() const { return y; }
/// @brief String representation of the object
/// @param verbose Verbosity level
@@ -102,7 +104,7 @@ public:
if (verbose > 0) {
msg << setw(14) << setfill(' ') << "dx [cm]" << ' ';
msg << setw(14) << setfill(' ') << "dy [cm]" << ' ';
- msg << setw(14) << setfill(' ') << "dxy [cm]" << ' ';
+ msg << setw(14) << setfill(' ') << "dxy [cm2]" << ' ';
msg << setw(14) << setfill(' ') << "dt [ns]" << ' ';
}
}
@@ -129,6 +131,7 @@ public:
int ExtIndex; ///< index of hit in the external branch
int IntIndex; ///< index of hit in the internal array
int iStation; ///< index of station in active stations array
+ int Det; ///< detector subsystem ID
double x; ///< x coordinate of position [cm]
double y; ///< y coordinate of position [cm]
double time; ///< hit time [ns]
@@ -136,8 +139,7 @@ public:
double dy; ///< y coordinate error [cm]
double dt; ///< time error [ns]
double dxy; ///< covariance between x and y [cm2]
- int Det; ///< detector subsystem ID
- L1Vector<int> mcPointIds {"CbmL1HitDebugInfo::mcPointIds"}; ///< indices of CbmL1MCPoint in L1->fvMCPoints array
+ int fPointID = -1; ///< index of matched MC point
};
#endif
diff --git a/reco/L1/CbmL1Performance.cxx b/reco/L1/CbmL1Performance.cxx
index 38dca67a4bacdfd37d772efbd6a43d7e4cad051e..2fce9259f9db132b1f6dfc5e9b082d3f84776c86 100644
--- a/reco/L1/CbmL1Performance.cxx
+++ b/reco/L1/CbmL1Performance.cxx
@@ -101,27 +101,23 @@ void CbmL1::TrackMatch()
map<int, int>& hitmap = prtra->hitMap; // how many hits from each mcTrack belong to current recoTrack
map<int, int> hitmapChain;
for (vector<int>::iterator ih = (prtra->Hits).begin(); ih != (prtra->Hits).end(); ++ih) {
-
- const int nMCPoints = fvHitDebugInfo[*ih].mcPointIds.size();
- for (int iP = 0; iP < nMCPoints; iP++) {
- int iMC = fvHitDebugInfo[*ih].mcPointIds[iP];
-
- // cout<<iMC<<" iMC"<<endl;
- int ID = -1;
- int chainID = -1;
- if (iMC >= 0) {
- ID = fvMCPoints[iMC].ID;
- chainID = fvMCTracks[ID].chainID;
- }
- if (hitmap.find(ID) == hitmap.end()) hitmap[ID] = 1;
- else {
- hitmap[ID] += 1;
- }
- if (hitmapChain.find(chainID) == hitmapChain.end()) hitmapChain[chainID] = 1;
- else {
- hitmapChain[chainID] += 1;
- }
- } // for iPoint
+ int iMC = fvHitDebugInfo[*ih].GetMCPointIndex();
+
+ // cout<<iMC<<" iMC"<<endl;
+ int ID = -1;
+ int chainID = -1;
+ if (iMC >= 0) {
+ ID = fvMCPoints[iMC].ID;
+ chainID = fvMCTracks[ID].chainID;
+ }
+ if (hitmap.find(ID) == hitmap.end()) hitmap[ID] = 1;
+ else {
+ hitmap[ID] += 1;
+ }
+ if (hitmapChain.find(chainID) == hitmapChain.end()) hitmapChain[chainID] = 1;
+ else {
+ hitmapChain[chainID] += 1;
+ }
} // for iHit
// RTrack <-> MCTrack identification
diff --git a/reco/L1/CbmL1ReadEvent.cxx b/reco/L1/CbmL1ReadEvent.cxx
index a87bd4748b8bd3874637031be5310a1013bea397..400efdeb0e6403826da9967973597253c1c82b25 100644
--- a/reco/L1/CbmL1ReadEvent.cxx
+++ b/reco/L1/CbmL1ReadEvent.cxx
@@ -1145,9 +1145,7 @@ void CbmL1::ReadEvent(CbmEvent* event)
fvExternalHits.reserve(nHits);
fvHitDebugInfo.reserve(nHits);
fvHitPointIndexes.reserve(nHits);
-
- fpIODataManager->ResetInputData();
- fpIODataManager->ReserveNhits(nHits);
+ fpIODataManager->ResetInputData(nHits);
fpIODataManager->SetNhitKeys(NStrips);
// ----- Fill
@@ -1461,13 +1459,8 @@ bool CbmL1::ReadMCPoint(CbmL1MCPoint* MC, int iPoint, int file, int event, int i
void CbmL1::HitMatch()
{
// Clear contents
- for (auto& hit : fvHitDebugInfo) {
- hit.mcPointIds.clear();
- }
-
- for (auto& point : fvMCPoints) {
- point.hitIds.clear();
- }
+ std::for_each(fvHitDebugInfo.begin(), fvHitDebugInfo.end(), [&](auto& hit) { hit.SetMCPointIndex(-1); });
+ std::for_each(fvMCPoints.begin(), fvMCPoints.end(), [&](auto& point) { point.hitIds.clear(); });
// Fill new contents
const int NHits = fvHitDebugInfo.size();
@@ -1475,7 +1468,7 @@ void CbmL1::HitMatch()
CbmL1HitDebugInfo& hit = fvHitDebugInfo[iH];
int iP = fvHitPointIndexes[iH];
if (iP >= 0) {
- hit.mcPointIds.push_back_no_warning(iP);
+ hit.SetMCPointIndex(iP);
fvMCPoints[iP].hitIds.push_back_no_warning(iH);
}
}
diff --git a/reco/L1/L1Algo/L1IODataManager.cxx b/reco/L1/L1Algo/L1IODataManager.cxx
index a76fad40bbde38b9338e4ff2882177a2e1619521..37913828d6ff099049c8188bfbac8c61c09d75cc 100644
--- a/reco/L1/L1Algo/L1IODataManager.cxx
+++ b/reco/L1/L1Algo/L1IODataManager.cxx
@@ -75,13 +75,14 @@ void L1IODataManager::ReadInputData(const std::string& fileName)
// ---------------------------------------------------------------------------------------------------------------------
//
-void L1IODataManager::ResetInputData() noexcept
+void L1IODataManager::ResetInputData(L1HitIndex_t nHits) noexcept
{
L1InputData tmp;
fInputData.Swap(tmp);
fLastStreamId = -1;
- fInputData.fStreamStartIndices.reserve(2000);
+ fInputData.fStreamStartIndices.reserve(2000); // TODO: What are these numbers? Please, put them into L1Constants.h
fInputData.fStreamStopIndices.reserve(2000);
+ fInputData.fHits.reserve(nHits);
}
// ---------------------------------------------------------------------------------------------------------------------
diff --git a/reco/L1/L1Algo/L1IODataManager.h b/reco/L1/L1Algo/L1IODataManager.h
index 49bcf6fa58008247f15ede4034cc394748f17e57..b7902481d6b81077285cf27457398ffd47f466de 100644
--- a/reco/L1/L1Algo/L1IODataManager.h
+++ b/reco/L1/L1Algo/L1IODataManager.h
@@ -57,8 +57,9 @@ public:
/// \note If one does not call this method, the underlying vector of hits will be filled with the time penalty
void ReserveNhits(L1HitIndex_t nHits) { fInputData.fHits.reserve(nHits); }
- /// Resets the input data block
- void ResetInputData() noexcept;
+ /// @brief Resets the input data block
+ /// @param nHits Number of hits to reserve
+ void ResetInputData(L1HitIndex_t nHits = 0) noexcept;
/// Pushes back a hit
/// \param hit An L1Hit object
diff --git a/reco/L1/L1Algo/L1InitManager.h b/reco/L1/L1Algo/L1InitManager.h
index 9a2f76cae719c6b4aa6e6d2de7e740d677a98da0..e85cc484ab3d252f128e7c6f6e659ec9c65ccafe 100644
--- a/reco/L1/L1Algo/L1InitManager.h
+++ b/reco/L1/L1Algo/L1InitManager.h
@@ -225,6 +225,11 @@ public:
/// \param filename Name of the output L1 parameters configuration
void SetOutputConfigName(const std::string& filename) { fConfigOutputName = filename; }
+ /// @brief Sets pseudo-random numbers generator seed
+ /// @param seed Seed value
+ /// @note The default seed is 0
+ void SetRandomSeed(unsigned int seed) { fParameters.fRandomSeed = seed; }
+
/// Sets target position
/// \param x Position X component [cm]
/// \param y Position Y component [cm]
diff --git a/reco/L1/L1Algo/L1Parameters.cxx b/reco/L1/L1Algo/L1Parameters.cxx
index 443b6660e6765e7c6c14dc2e322949a0de87e254..396c09b151a59ec9466bb4defa63510a353b853b 100644
--- a/reco/L1/L1Algo/L1Parameters.cxx
+++ b/reco/L1/L1Algo/L1Parameters.cxx
@@ -232,6 +232,7 @@ std::string L1Parameters::ToString(int verbosity, int indentLevel) const
aStream << indent << indentChar << "Max number of threads: " << L1Constants::size::kMaxNthreads << '\n';
aStream << indent << indentChar << "Max number of triplets: " << L1Constants::size::kMaxNtriplets << '\n';
aStream << indent << kCLb << "RUNTIME CONSTANTS:\n" << kCL;
+ aStream << indent << indentChar << "Random seed: " << fRandomSeed << '\n';
aStream << indent << indentChar << "Max number of doublets per singlet: " << fMaxDoubletsPerSinglet << '\n';
aStream << indent << indentChar << "Max number of triplets per doublet: " << fMaxTripletPerDoublets << '\n';
aStream << indent << indentChar << "Tracking level: " << fTrackingLevel << '\n';
diff --git a/reco/L1/L1Algo/L1Parameters.h b/reco/L1/L1Algo/L1Parameters.h
index bfe7ba29142d5da28b9a899771d436a77c8b8f7d..03c235f250053cd887c93cc1800ba6f1932fb5ab 100644
--- a/reco/L1/L1Algo/L1Parameters.h
+++ b/reco/L1/L1Algo/L1Parameters.h
@@ -71,12 +71,21 @@ public:
/// String representation of the class contents
std::string ToString(int verbosity = 0, int indentLevel = 0) const;
+ // TODO: SZh 15.05.2023: Setters should be removed from the class, since parameters
+ // can be initialized only via L1InitManager.
/// Sets upper-bound cut on max number of doublets per one singlet
void SetMaxDoubletsPerSinglet(unsigned int value) { fMaxDoubletsPerSinglet = value; }
/// Sets upper-bound cut on max number of triplets per one doublet
void SetMaxTripletPerDoublets(unsigned int value) { fMaxTripletPerDoublets = value; }
+ /// @brief Gets first active station index within detector subsystem
+ /// @param detectorID Detector ID
+ int GetFirstActiveStationIndex(L1DetectorID detectorID) const
+ {
+ return std::accumulate(fNstationsActive.cbegin(), fNstationsActive.cbegin() + static_cast<int>(detectorID), 0);
+ }
+
/// Gets upper-bound cut on max number of doublets per one singlet
unsigned int GetMaxDoubletsPerSinglet() const { return fMaxDoubletsPerSinglet; }
@@ -181,6 +190,12 @@ public:
/// Gets L1FieldValue object at primary vertex
const L1FieldValue& GetVertexFieldValue() const { return fVertexFieldValue; }
+ /// @brief Gets random seed
+ ///
+ /// If random seed is zero, std::random_device is used to seed the random number generator.
+ int GetRandomSeed() const { return fRandomSeed; }
+
+
/// Gets ghost suppression flag
int GetGhostSuppression() const { return fGhostSuppression; }
@@ -268,6 +283,8 @@ private:
int fTrackingLevel = 0; ///< tracking level
int fGhostSuppression = 0; ///< flag: if true, ghost tracks are suppressed
float fMomentumCutOff = 0; ///< minimum momentum of tracks
+ int fRandomSeed = 1; ///< random seed
+
// ***************************
// ** Flags for development **
@@ -278,7 +295,7 @@ private:
bool fDevIsMatchDoubletsViaMc {false}; ///< Flag to match doublets using MC information
bool fDevIsMatchTripletsViaMc {false}; ///< Flag to match triplets using Mc information
bool fDevIsExtendTracksViaMc {false}; ///< Flag to extend tracks using Mc information
- bool fDevIsSuppressOverlapHitsViaMc {false}; ///< Flag to match hits in overlaps using Mc information
+ bool fDevIsSuppressOverlapHitsViaMc {false}; ///< Flag to match hits in overlaps using Mc information
bool fDevIsParSearchWUsed = false; ///< Flag: when true, the parametrized search windows are used in track
///< finder; when false, the Kalman filter is used instead
@@ -310,6 +327,7 @@ private:
ar& fTrackingLevel;
ar& fGhostSuppression;
ar& fMomentumCutOff;
+ ar& fRandomSeed;
ar& fDevIsIgnoreHitSearchAreas;
ar& fDevIsUseOfOriginalField;
diff --git a/reco/L1/L1Algo/L1Station.h b/reco/L1/L1Algo/L1Station.h
index 0b86797812e65485c35e5981bab2afaff63d73dd..f22946a8975fed64034e2c43bfe355daecc3c390 100644
--- a/reco/L1/L1Algo/L1Station.h
+++ b/reco/L1/L1Algo/L1Station.h
@@ -21,6 +21,9 @@
///
class L1Station {
public:
+ // TODO: SZh 12.05.2022: Rewrite type into L1DetectorID, change detector indexing scheme
+ // TODO: SZh 12.05.2022: Provide getters to stations
+
int type = undef::kI32;
int timeInfo = undef::kI32; ///< flag: if time information can be used
int fieldStatus = undef::kI32; ///< flag: 1 - station is INSIDE the field, 0 - station is OUTSIDE the field
diff --git a/reco/L1/L1Algo/L1Utils.h b/reco/L1/L1Algo/L1Utils.h
index bd1b798ab9c0545d14901a3520566f86b9f36bc8..41ad13f8f5ee56ad8ac54a610224beb72dc94f3f 100644
--- a/reco/L1/L1Algo/L1Utils.h
+++ b/reco/L1/L1Algo/L1Utils.h
@@ -26,11 +26,14 @@
#include "CaSimd.h"
#include "L1Def.h"
+// TODO: SZh 16.05.2023: Split this class into several headers and place them into L1Algo::utils
+
/// Class contains some utility functions for L1Algo
namespace L1Utils
{
/// NaN value for float
+ // TODO: SZh 16.05.2023: Replace this kNaN with L1Undef::kF
constexpr float kNaN {std::numeric_limits<float>::signaling_NaN()};
/// Comparison method for floats
diff --git a/reco/L1/L1Algo/utils/CaAlgoRandom.cxx b/reco/L1/L1Algo/utils/CaAlgoRandom.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..218abc50bcf55c1f94163b24d98839360bb867a4
--- /dev/null
+++ b/reco/L1/L1Algo/utils/CaAlgoRandom.cxx
@@ -0,0 +1,34 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// @file CaAlgoRandom.cxx
+/// @brief Random generator utility class for CA tracking (implementation)
+/// @since 16.05.2023
+/// @author Sergei Zharko <s.zharko@gsi.de>
+
+#include "CaAlgoRandom.h"
+
+using ca::algo::Random;
+
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+Random::Random() : Random(1) {}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+Random::Random(int seed) { this->SetSeed(seed); }
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void Random::SetSeed(int seed)
+{
+ if (seed) { fSeed = seed; }
+ else {
+ std::random_device rd;
+ unsigned int randomDeviceSeed = rd();
+ fSeed = randomDeviceSeed;
+ }
+ fGenerator.seed(fSeed);
+}
diff --git a/reco/L1/L1Algo/utils/CaAlgoRandom.h b/reco/L1/L1Algo/utils/CaAlgoRandom.h
new file mode 100644
index 0000000000000000000000000000000000000000..a98d18e850e5f3a4237e1a4988807ba711561ebf
--- /dev/null
+++ b/reco/L1/L1Algo/utils/CaAlgoRandom.h
@@ -0,0 +1,94 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// @file CaAlgoRandom.h
+/// @brief Random generator utility class for CA tracking
+/// @since 16.05.2023
+/// @authors Sergei Zharko <s.zharko@gsi.de>, Sergey Gorbunov <se.gorbunov@gsi.de>
+
+#ifndef CaAlgoRandom_h
+#define CaAlgoRandom_h 1
+
+#include <cstdint>
+#include <random>
+#include <type_traits>
+
+namespace ca::algo
+{
+ /// @brief A class, providing ROOT-free access to randomly generated variables
+ class Random {
+ using GeneratorType_t = std::mt19937_64;
+
+ public:
+ /// @brief Default constructor
+ ///
+ /// Sets seed equal to 1
+ Random();
+
+ /// @brief Constructor
+ /// @param seed A random seed
+ ///
+ /// If 0 is selected for the random seed, the generator is seeded with std::random_device
+ Random(int seed);
+
+ /// @brief Destructor
+ ~Random() = default;
+
+ /// @brief Copy constructor
+ Random(const Random& other) = delete;
+
+ /// @brief Move constructor
+ Random(Random&& other) = delete;
+
+ /// @brief Copy assignment operator
+ Random& operator=(const Random&) = delete;
+
+ /// @brief Move assignment operator
+ Random& operator=(Random&&) = delete;
+
+ /// @brief Gets seed of the random generator
+ int GetSeed() const { return fSeed; }
+
+ /// @brief Sets seed to the random generator
+ /// @param seed A random seed
+ ///
+ /// If 0 is selected for the random seed, the generator is seeded with std::random_device
+ void SetSeed(int seed);
+
+ /// @brief Returns a normally distributed random value, limited within a selected sigma range
+ /// @tparam T Type of floating point
+ /// @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;
+
+ private:
+ mutable GeneratorType_t fGenerator; ///< Random number generator
+ unsigned int fSeed = 1; ///< Random number seed
+ };
+} // namespace ca::algo
+
+
+// *********************************************************
+// ** Template and inline function implementation **
+// *********************************************************
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+template<typename T, std::enable_if_t<std::is_floating_point<T>::value, bool>>
+T ca::algo::Random::BoundedGaus(const T& mean, const T& sigma, const T& nSigmas) const
+{
+ assert(sigma > 0 && std::isfinite(sigma));
+ assert(nSigmas > 0 && std::isfinite(nSigmas));
+ std::normal_distribution rndGaus {mean, sigma};
+ double res = 0;
+ do {
+ res = rndGaus(fGenerator);
+ } while (std::fabs(res - mean) > nSigmas * sigma);
+ return res;
+}
+
+
+#endif // CaAlgoRandom_h
diff --git a/reco/L1/L1Algo/utils/L1AlgoEfficiencyPerformance.h b/reco/L1/L1Algo/utils/L1AlgoEfficiencyPerformance.h
index 512d6ce217825b0178ccc4f9148844f67dc288d7..6d767fe1040ac1c4b465a0f4bd879c4af03c6018 100644
--- a/reco/L1/L1Algo/utils/L1AlgoEfficiencyPerformance.h
+++ b/reco/L1/L1Algo/utils/L1AlgoEfficiencyPerformance.h
@@ -224,12 +224,11 @@ bool L1AlgoEfficiencyPerformance<NHits>::AddOne(L1HitIndex_t* iHits)
// obtain mc data
for (int iih = 0; iih < NHits; iih++) {
- int nMC = fL1->fvHitDebugInfo[iHits[iih]].mcPointIds.size();
- for (int iMC = 0; iMC < nMC; iMC++) {
- const int iMCP = fL1->fvHitDebugInfo[iHits[iih]].mcPointIds[iMC];
+ const int iMCP = fL1->fvHitDebugInfo[iHits[iih]].GetMCPointIndex();
+ if (iMCP > -1) {
int mcId = fL1->fvMCPoints[iMCP].ID;
mcIds[iih].push_back(mcId);
- } // for mcPoints
+ }
} // for hits
// find all reconstructed track id-s
@@ -250,11 +249,14 @@ bool L1AlgoEfficiencyPerformance<NHits>::AddOne(L1HitIndex_t* iHits)
}
// use first found // TODO: save all!!!
+ // UPD: SZh 24.05.2023: We save only one point index (the best matched)
vector<int>& mcsN = mcIds[NHits - 1];
for (unsigned int i = 0; i < mcsN.size(); i++) {
if (mcsN[i] >= 0) {
trlet.mcTrackId = mcsN[i];
- trlet.iStation = fL1->fvMCPoints[fL1->fvHitDebugInfo[iHits[0]].mcPointIds[0]].iStation;
+ int iMCP = fL1->fvHitDebugInfo[iHits[0]].GetMCPointIndex();
+ assert(iMCP > -1);
+ trlet.iStation = fL1->fvMCPoints[fL1->fvHitDebugInfo[iHits[0]].GetMCPointIndex()].iStation;
break;
}
}
diff --git a/reco/L1/L1LinkDef.h b/reco/L1/L1LinkDef.h
index ad9b390b8c78b25641fe023de28ebe4793ff7447..e5d1a2a9478c20f88431791c829bf07ef12bcfc2 100644
--- a/reco/L1/L1LinkDef.h
+++ b/reco/L1/L1LinkDef.h
@@ -35,7 +35,7 @@
#pragma link C++ class CbmCaInputQaTrd + ;
#pragma link C++ class CbmCaInputQaTof + ;
#pragma link C++ enum cbm::ca::ETrackType;
-#pragma link C++ enum L1DetectorID;
+#pragma link C++ enum class L1DetectorID;
#pragma link C++ class cbm::ca::OutputQa + ;
#pragma link C++ class ca::tools::WindowFinder + ;
diff --git a/reco/L1/catools/CaToolsHitRecord.h b/reco/L1/catools/CaToolsHitRecord.h
new file mode 100644
index 0000000000000000000000000000000000000000..9915a0661b8c7cef1330971d8fc3cee0de335b17
--- /dev/null
+++ b/reco/L1/catools/CaToolsHitRecord.h
@@ -0,0 +1,45 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// \file CbmCaHitRecord.h
+/// \brief A header for cbm::ca::HitRecord structure
+/// \since 15.05.2023
+/// \author S.Zharko <s.zharko@gsi.de>
+
+#ifndef CbmCaHitRecord_h
+#define CbmCaHitRecord_h 1
+
+#include <cstdint>
+
+namespace ca::tools
+{
+ /// @brief A helper structure to store hits information from different detectors in a uniform manner
+ ///
+ /// @note This structure is to be used ONLY at the TS/event initialization stage to fill hit arrays of different
+ /// types.
+ ///
+ struct HitRecord {
+ double fX = -1.; ///< x component of hit position [cm]
+ double fY = -1.; ///< y component of hit position [cm]
+ double fDx = -1.; ///< error of x component of hit position [cm]
+ double fDy = -1.; ///< error of y component of hit position [cm]
+ double fDxy = -1.; ///< correlation between x and y components [cm]
+ double fU = -1.; ///< hit position in direction of front strips [cm]
+ double fV = -1.; ///< hit position in direction of back strips [cm]
+ double fDu = -1.; ///< hit position error in direction of front strips [cm]
+ double fDv = -1.; ///< hit position error in direction of back strips [cm]
+ double fZ = -1.; ///< z component of hit position [cm]
+ double fT = -1.; ///< time of hit [ns]
+ double fDt = -1.; ///< time error of hit [ns]
+ int64_t fDataStream = -1; ///< Global index of detector module
+ int fPointId = -1; ///< index of MC point
+ int fExtId = -2; ///< external index of hit
+ int fStaId = -2; ///< index of active tracking station
+ int fStripF = -2; ///< index of front strip
+ int fStripB = -2; ///< index of back strip
+ int fDet = -2; ///< detector ID
+ };
+} // namespace ca::tools
+
+#endif // CbmCaHitRecord_h
diff --git a/reco/L1/catools/CaToolsMCData.h b/reco/L1/catools/CaToolsMCData.h
index 72c3d95e07ff5d39f9812df18822839dd0c03678..a20946c73b5a94dc8fc05cc5feb10b9622681b54 100644
--- a/reco/L1/catools/CaToolsMCData.h
+++ b/reco/L1/catools/CaToolsMCData.h
@@ -96,15 +96,15 @@ namespace ca::tools
return std::accumulate(fvNofPointsUsed.cbegin(), fvNofPointsUsed.cbegin() + static_cast<int>(detID), 0);
}
- /// @brief Gets the first point index for a given detector subsystem
+ /// @brief Gets the next index of point, which is expected being after the last one for a given detector
/// @param detID Detector ID
int GetLastPointIndex(L1DetectorID detID) const
{
- return std::accumulate(fvNofPointsUsed.cbegin(), fvNofPointsUsed.cbegin() + static_cast<int>(detID) + 1, 0) - 1;
+ return std::accumulate(fvNofPointsUsed.cbegin(), fvNofPointsUsed.cbegin() + static_cast<int>(detID) + 1, 0);
}
/// @brief Calculates global index of MC point
- /// @param iPointLocal Local index of MC point
+ /// @param iPointLocal Local index of MC psinoint
/// @param detID Detector ID
///
/// The function calculates global external index of MC point as a sum of a given local index and total provided
@@ -114,7 +114,6 @@ namespace ca::tools
return iPointLocal + std::accumulate(fvNofPointsOrig.cbegin(), fvNofPointsOrig.cbegin() + int(detID), 0);
}
-
/// Gets number of tracks in this event/TS
int GetNofTracks() const { return fvTracks.size(); }
diff --git a/reco/L1/catools/CaToolsMCPoint.h b/reco/L1/catools/CaToolsMCPoint.h
index 158f35e6e3fb35e7377a718f94ca9bdfb790968e..5b5466dcd2398e78a13079a3e222bb3553b05ef9 100644
--- a/reco/L1/catools/CaToolsMCPoint.h
+++ b/reco/L1/catools/CaToolsMCPoint.h
@@ -62,6 +62,11 @@ namespace ca::tools
/// @brief Gets MC event ID
int GetEventId() const { return fLinkKey.fEvent; }
+ /// @brief Gets MC external ID
+ /// @note External ID of a point shifted from the original external ID by the total number of points,
+ /// produced in the preceding detector subsystems.
+ int GetExternalId() const { return fLinkKey.fIndex; }
+
/// @brief Gets MC file ID
int GetFileId() const { return fLinkKey.fFile; }
@@ -350,6 +355,7 @@ namespace ca::tools
L1DetectorID fDetectorId; ///< Detector ID of MC point
+ // TODO: SZh 17.05.2023: Check, if there are more then one index can be added
L1Vector<int> fvHitIndexes {"ca::tools::MCPoint::fvHitIndexes"}; ///< Indexes of hits, assigned to this point
};
} // namespace ca::tools
diff --git a/reco/L1/qa/CbmCaOutputQa.cxx b/reco/L1/qa/CbmCaOutputQa.cxx
index 69f1a33db568823fb803cedd34b75360279eb8d0..f2dc7dab18a7b7c997fdd2326aff5c1e8f9a300f 100644
--- a/reco/L1/qa/CbmCaOutputQa.cxx
+++ b/reco/L1/qa/CbmCaOutputQa.cxx
@@ -29,6 +29,11 @@ using cbm::ca::OutputQa;
//
OutputQa::OutputQa(int verbose, bool isMCUsed) : CbmQaTask("CbmCaOutputQa", "caout", verbose, isMCUsed)
{
+ // Create MC module
+
+ // Create TS reader
+ fpTSReader = std::make_unique<TimeSliceReader>();
+
// Turn on default track classes
AddTrackType(ETrackType::kAll);
AddTrackType(ETrackType::kGhost);
@@ -451,7 +456,7 @@ InitStatus OutputQa::InitDataBranches()
if (!fpDataManager.get()) { fpDataManager = std::make_shared<L1IODataManager>(); }
// Initialize time slice reader instance
- if (!fpTSReader.get()) { fpTSReader = std::make_unique<TimeSliceReader>(fTrackingMode); }
+ fpTSReader->SetTrackingMode(fTrackingMode);
fpTSReader->SetDetector(L1DetectorID::kMvd, fbUseMvd);
fpTSReader->SetDetector(L1DetectorID::kSts, fbUseSts);
fpTSReader->SetDetector(L1DetectorID::kMuch, fbUseMuch);
@@ -462,12 +467,11 @@ InitStatus OutputQa::InitDataBranches()
fpTSReader->RegisterTracksContainer(fvRecoTracks);
fpTSReader->RegisterQaHitContainer(fvHits);
fpTSReader->RegisterHitIndexContainer(fvHitIds);
-
if (!fpTSReader->InitRun()) { return kFATAL; }
// Initialize MC module
if (IsMCUsed()) {
- if (!fpMCModule.get()) { fpMCModule = std::make_unique<MCModule>(fVerbose, fPerformanceMode); }
+ fpMCModule = std::make_shared<MCModule>(fVerbose, fPerformanceMode);
fpMCModule->SetDetector(L1DetectorID::kMvd, fbUseMvd);
fpMCModule->SetDetector(L1DetectorID::kSts, fbUseSts);
fpMCModule->SetDetector(L1DetectorID::kMuch, fbUseMuch);
@@ -480,9 +484,9 @@ InitStatus OutputQa::InitDataBranches()
fpMCModule->RegisterQaHitContainer(fvHits);
fpMCModule->RegisterParameters(fpParameters);
fpMCModule->RegisterFirstHitIndexes(fpTSReader->GetHitFirstIndexDet());
-
if (!fpMCModule->InitRun()) { return kFATAL; }
}
+
return kSUCCESS;
}
@@ -730,38 +734,20 @@ InitStatus OutputQa::InitTimeSlice()
int nHits = 0;
int nRecoTracks = 0;
- // Read MC tracks and points
- if (IsMCUsed()) {
- fpMCModule->InitEvent(nullptr);
- nMCPoints = fMCData.GetNofPoints();
- nMCTracks = fMCData.GetNofTracks();
- }
// Read reconstructed input
fpTSReader->InitTimeSlice();
nHits = fvHits.size();
nRecoTracks = fvRecoTracks.size();
- static bool bDo = true;
- if (bDo) {
- for (int iP = 0; iP < nMCPoints; ++iP) {
- const auto& point = fMCData.GetPoint(iP);
- LOG(info) << iP << ' ' << (int) point.GetDetectorId() << ' ' << point.GetStationId();
- }
-
- for (int iH = 0; iH < nHits; ++iH) {
- const auto& hit = fvHits[iH];
- LOG(info) << iH << ' ' << (int) hit.GetDetectorType() << ' ' << hit.GetStationId();
- }
- bDo = false;
- }
-
- for (const auto& hit : fvHits) {
- assert(hit.GetMCPointIndexes().size() < 2);
- }
-
// Match tracks and points
- if (IsMCUsed()) { fpMCModule->MatchRecoAndMC(); }
+ // Read MC tracks and points
+ if (IsMCUsed()) {
+ fpMCModule->InitEvent(nullptr);
+ nMCPoints = fMCData.GetNofPoints();
+ nMCTracks = fMCData.GetNofTracks();
+ fpMCModule->MatchRecoAndMC();
+ }
LOG_IF(info, fVerbose > 1) << fName << ": Data sample consists of " << nHits << " hits, " << nRecoTracks
<< " reco tracks, " << nMCTracks << " MC tracks, " << nMCPoints << " MC points";
diff --git a/reco/L1/qa/CbmCaOutputQa.h b/reco/L1/qa/CbmCaOutputQa.h
index 1cc6216224126eb41189fe924c0dd1b9e9bfdea6..e3902ae4ff1e1564cbfd32a4f48c04043c5988bf 100644
--- a/reco/L1/qa/CbmCaOutputQa.h
+++ b/reco/L1/qa/CbmCaOutputQa.h
@@ -154,10 +154,10 @@ namespace cbm::ca
void SetUseTof(bool flag = true) { fbUseTof = flag; }
/// @brief Sets STS tracking mode
- void SetStsTrackingMode() { fTrackingMode = ECbmTrackingMode::kSTS; }
+ void SetStsTrackingMode() { fTrackingMode = ECbmCaTrackingMode::kSTS; }
/// @brief Sets mCBM global tracking mode
- void SetMcbmTrackingMode() { fTrackingMode = ECbmTrackingMode::kMCBM; }
+ void SetMcbmTrackingMode() { fTrackingMode = ECbmCaTrackingMode::kMCBM; }
/// @brief Sets performance mode
/// @param mode Performance mode (3 is default: TODO - test)
@@ -232,10 +232,10 @@ namespace cbm::ca
int fPerformanceMode = 3; ///< Performance mode
- ECbmTrackingMode fTrackingMode = ECbmTrackingMode::kSTS; ///< Tracking mode
+ ECbmCaTrackingMode fTrackingMode = ECbmCaTrackingMode::kSTS; ///< Tracking mode
std::unique_ptr<TimeSliceReader> fpTSReader = nullptr; ///< Reader of the time slice
- std::unique_ptr<cbm::ca::MCModule> fpMCModule = nullptr; ///< MC module
+ std::shared_ptr<MCModule> fpMCModule = nullptr; ///< MC module
std::shared_ptr<L1IODataManager> fpDataManager = nullptr; ///< Data manager
std::shared_ptr<::ca::tools::Debugger> fpDebugger = nullptr; ///< Debugger
std::shared_ptr<L1Parameters> fpParameters = nullptr; ///< Tracking parameters object
diff --git a/reco/L1/qa/CbmCaTrackFitQa.h b/reco/L1/qa/CbmCaTrackFitQa.h
index 91686b14eabac1a826719f4fbefd5db0913ede19..051d200ba59d30d7ad47355d48c18fa62ce64666 100644
--- a/reco/L1/qa/CbmCaTrackFitQa.h
+++ b/reco/L1/qa/CbmCaTrackFitQa.h
@@ -183,19 +183,19 @@ namespace cbm::ca
/// @param weight Weight
void FillResAndPull(ETrackParType type, double recoVal, double recoErr, double trueVal, double w);
- TrackParArray_t<TH1F*> fvphResiduals = {0}; ///< Residuals for different track parameters
- TrackParArray_t<TH1F*> fvphPulls = {0}; ///< Pulls for different track parameters
+ TrackParArray_t<TH1F*> fvphResiduals = {{0}}; ///< Residuals for different track parameters
+ TrackParArray_t<TH1F*> fvphPulls = {{0}}; ///< Pulls for different track parameters
- TrackParArray_t<bool> fvbParIgnored = {0}; ///< Flag: true - parameter is ignored
+ TrackParArray_t<bool> fvbParIgnored = {{0}}; ///< Flag: true - parameter is ignored
// ** Distribution properties **
- TrackParArray_t<int> fvRBins = {0}; ///< Number of bins, residuals
- TrackParArray_t<double> fvRLo = {0}; ///< Lower boundary, residuals
- TrackParArray_t<double> fvRUp = {0}; ///< Upper boundary, residuals
+ TrackParArray_t<int> fvRBins = {{0}}; ///< Number of bins, residuals
+ TrackParArray_t<double> fvRLo = {{0}}; ///< Lower boundary, residuals
+ TrackParArray_t<double> fvRUp = {{0}}; ///< Upper boundary, residuals
- TrackParArray_t<int> fvPBins = {0}; ///< Number of bins, pulls
- TrackParArray_t<double> fvPLo = {0}; ///< Lower boundary, pulls
- TrackParArray_t<double> fvPUp = {0}; ///< Upper boundary, pulls
+ TrackParArray_t<int> fvPBins = {{0}}; ///< Number of bins, pulls
+ TrackParArray_t<double> fvPLo = {{0}}; ///< Lower boundary, pulls
+ TrackParArray_t<double> fvPUp = {{0}}; ///< Upper boundary, pulls
TString fsTitle = ""; ///< Title of the point