diff --git a/core/detectors/much/CbmMuchTrackingInterface.h b/core/detectors/much/CbmMuchTrackingInterface.h
index 5ba7677dba529eaf0de04b4a55303222a0f86ac5..ec82756563560582aa386bb0535e5063b561031e 100644
--- a/core/detectors/much/CbmMuchTrackingInterface.h
+++ b/core/detectors/much/CbmMuchTrackingInterface.h
@@ -95,10 +95,7 @@ public:
/// Gets a tracking station by the address of element
/// \param address Unique element address
/// \return Local index of the tracking station
- int GetTrackingStationIndex(int address) const
- {
- return CbmMuchGeoScheme::GetStationIndex(address) * 3 + CbmMuchGeoScheme::GetLayerIndex(address);
- }
+ int GetTrackingStationIndex(int address) const { return CbmMuchGeoScheme::GetStationIndex(address); }
/// Gets max size of a station along the X-axis
/// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
diff --git a/macro/mcbm/mcbm_qa.C b/macro/mcbm/mcbm_qa.C
index 64336d6abf83b82a2bfe55cb486be1579b6b36a1..febe63f12e5918df6a35ced8748b561265dab5da 100644
--- a/macro/mcbm/mcbm_qa.C
+++ b/macro/mcbm/mcbm_qa.C
@@ -76,7 +76,7 @@ void mcbm_qa(Int_t nEvents = 0, TString dataset = "data/mcbm_beam_2020_03_test",
// ------------------------------------------------------------------------
// ----- Some global switches -----------------------------------------
- if (setupName == "mcbm_beam_2022_05_23_nickel") { setup->RemoveModule(ECbmModuleId::kMuch); }
+ //if (setupName == "mcbm_beam_2022_05_23_nickel") { setup->RemoveModule(ECbmModuleId::kMuch); }
//bool eventBased = !sEvBuildRaw.IsNull();
bool bUseMvd = setup->IsActive(ECbmModuleId::kMvd);
diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx
index 38376d99dbe74465c9886ab8fb8ff812bae567e1..4cf54dc1b3add62c736fe771400da8b0de133f83 100644
--- a/reco/L1/CbmL1.cxx
+++ b/reco/L1/CbmL1.cxx
@@ -362,7 +362,9 @@ InitStatus CbmL1::Init()
}
// *****************************
- // ** Geometry initialization **
+ // ** **
+ // ** GEOMETRY INITIALIZATION **
+ // ** **
// *****************************
// Read parameters object from a binary file
@@ -376,6 +378,8 @@ InitStatus CbmL1::Init()
fNTofStationsGeom = fInitManager.GetNstationsGeometry(L1DetectorID::kTof);
fNStationsGeom = fInitManager.GetNstationsGeometry();
+ // TODO: SZh 15.08.2023: This particular hack is related only to the mcbm_beam_2021_07_surveyed setup.
+ // Shell we still support it?
if (fMissingHits) { CbmTofTrackingInterface::Instance()->FixHitsStationsMismatch(); }
}
// Parameters initialization, based on the FairRunAna chain
@@ -448,22 +452,9 @@ InitStatus CbmL1::Init()
if (fUseTOF) { vActiveTrackingDetectorIDs.insert(L1DetectorID::kTof); }
fInitManager.SetActiveDetectorIDs(vActiveTrackingDetectorIDs);
- // *********************************************************************
- // ** Counting numbers of stations for different detector subsystems **
- // *********************************************************************
-
-
- // Provide crosscheck number of stations for the fInitManagera
- fInitManager.SetNstations(L1DetectorID::kMvd, fNMvdStationsGeom);
- fInitManager.SetNstations(L1DetectorID::kSts, fNStsStationsGeom);
- fInitManager.SetNstations(L1DetectorID::kMuch, fNMuchStationsGeom);
- fInitManager.SetNstations(L1DetectorID::kTrd, fNTrdStationsGeom);
- fInitManager.SetNstations(L1DetectorID::kTof, fNTofStationsGeom);
-
-
- // ***************************************
- // ** Stations geometry initialization **
- // ***************************************
+ // *************************************
+ // ** Stations layout initialization **
+ // *************************************
std::vector<L1BaseStationInfo> vStations;
vStations.reserve(100);
@@ -472,6 +463,7 @@ InitStatus CbmL1::Init()
if (fUseMVD) {
for (int iSt = 0; iSt < fNMvdStationsGeom; ++iSt) {
auto stationInfo = L1BaseStationInfo(L1DetectorID::kMvd, iSt);
+ // TODO: SZh 15.08.2022: Replace station type with L1DetectorID
stationInfo.SetStationType(1); // MVD
stationInfo.SetTimeInfo(mvdInterface->IsTimeInfoProvided(iSt));
stationInfo.SetFieldStatus(fTrackingMode == L1Algo::TrackingMode::kMcbm ? 0 : 1);
@@ -488,8 +480,8 @@ InitStatus CbmL1::Init()
if (fvmDisabledStationIDs[L1DetectorID::kMvd].find(iSt) != fvmDisabledStationIDs[L1DetectorID::kMvd].cend()) {
stationInfo.SetTrackingStatus(false);
}
+ fInitManager.AddStation(stationInfo);
LOG(info) << "- MVD station " << iSt << " at z = " << stationInfo.GetZdouble() << " cm ";
- vStations.push_back(stationInfo);
}
}
@@ -497,6 +489,7 @@ InitStatus CbmL1::Init()
if (fUseSTS) {
for (int iSt = 0; iSt < fNStsStationsGeom; ++iSt) {
auto stationInfo = L1BaseStationInfo(L1DetectorID::kSts, iSt);
+ // TODO: SZh 15.08.2022: Replace station type with L1DetectorID
stationInfo.SetStationType(0); // STS
stationInfo.SetTimeInfo(stsInterface->IsTimeInfoProvided(iSt));
stationInfo.SetTimeInfo(L1Algo::TrackingMode::kMcbm != fTrackingMode ? stsInterface->IsTimeInfoProvided(iSt)
@@ -516,8 +509,8 @@ InitStatus CbmL1::Init()
if (fvmDisabledStationIDs[L1DetectorID::kSts].find(iSt) != fvmDisabledStationIDs[L1DetectorID::kSts].cend()) {
stationInfo.SetTrackingStatus(false);
}
+ fInitManager.AddStation(stationInfo);
LOG(info) << "- STS station " << iSt << " at z = " << stationInfo.GetZdouble() << " cm ";
- vStations.push_back(stationInfo);
}
}
@@ -525,6 +518,7 @@ InitStatus CbmL1::Init()
if (fUseMUCH) {
for (int iSt = 0; iSt < fNMuchStationsGeom; ++iSt) {
auto stationInfo = L1BaseStationInfo(L1DetectorID::kMuch, iSt);
+ // TODO: SZh 15.08.2022: Replace station type with L1DetectorID
stationInfo.SetStationType(2); // MuCh
stationInfo.SetTimeInfo(muchInterface->IsTimeInfoProvided(iSt));
stationInfo.SetTimeInfo(L1Algo::TrackingMode::kMcbm != fTrackingMode ? stsInterface->IsTimeInfoProvided(iSt)
@@ -543,7 +537,7 @@ InitStatus CbmL1::Init()
if (fvmDisabledStationIDs[L1DetectorID::kMuch].find(iSt) != fvmDisabledStationIDs[L1DetectorID::kMuch].cend()) {
stationInfo.SetTrackingStatus(false);
}
- vStations.push_back(stationInfo);
+ fInitManager.AddStation(stationInfo);
LOG(info) << "- MuCh station " << iSt << " at z = " << stationInfo.GetZdouble() << " cm";
}
}
@@ -552,6 +546,7 @@ InitStatus CbmL1::Init()
if (fUseTRD) {
for (int iSt = 0; iSt < fNTrdStationsGeom; ++iSt) {
auto stationInfo = L1BaseStationInfo(L1DetectorID::kTrd, iSt);
+ // TODO: SZh 15.08.2022: Replace station type with L1DetectorID
stationInfo.SetStationType((iSt == 1 || iSt == 3) ? 6 : 3); // MuCh
stationInfo.SetTimeInfo(trdInterface->IsTimeInfoProvided(iSt));
stationInfo.SetTimeInfo(L1Algo::TrackingMode::kMcbm != fTrackingMode ? stsInterface->IsTimeInfoProvided(iSt)
@@ -574,7 +569,7 @@ InitStatus CbmL1::Init()
if (fvmDisabledStationIDs[L1DetectorID::kTrd].find(iSt) != fvmDisabledStationIDs[L1DetectorID::kTrd].cend()) {
stationInfo.SetTrackingStatus(false);
}
- vStations.push_back(stationInfo);
+ fInitManager.AddStation(stationInfo);
LOG(info) << "- TRD station " << iSt << " at z = " << stationInfo.GetZdouble() << " cm";
}
}
@@ -583,6 +578,7 @@ InitStatus CbmL1::Init()
if (fUseTOF) {
for (int iSt = 0; iSt < fNTofStationsGeom; ++iSt) {
auto stationInfo = L1BaseStationInfo(L1DetectorID::kTof, iSt);
+ // TODO: SZh 15.08.2022: Replace station type with L1DetectorID
stationInfo.SetStationType(4);
stationInfo.SetTimeInfo(tofInterface->IsTimeInfoProvided(iSt));
stationInfo.SetTimeInfo(L1Algo::TrackingMode::kMcbm != fTrackingMode ? stsInterface->IsTimeInfoProvided(iSt)
@@ -602,74 +598,19 @@ InitStatus CbmL1::Init()
if (fvmDisabledStationIDs[L1DetectorID::kTof].find(iSt) != fvmDisabledStationIDs[L1DetectorID::kTof].cend()) {
stationInfo.SetTrackingStatus(false);
}
- vStations.push_back(stationInfo);
+ fInitManager.AddStation(stationInfo);
LOG(info) << "- TOF station " << iSt << " at z = " << stationInfo.GetZdouble() << " cm";
}
}
- // ** Order stations in Z
- {
- std::sort(vStations.begin(), vStations.end());
- }
+ // Init station layout: sort stations in z-position and init maps of station local, geo and active indices
+ fInitManager.InitStationLayout();
// ****************************************
- // ** **
- // ** Initialize material maps **
- // ** **
+ // ** Material maps initialization **
// ****************************************
+ this->GenerateMaterialMaps();
- {
- LOG(info) << "- Generate material maps..";
- auto timerStart = std::chrono::high_resolution_clock::now();
-
- ca::tools::MaterialHelper matHelper;
- matHelper.SetSafeMaterialInitialization(fDoSafeMaterialInitialization);
-
- if (!fMatBudgetParallelProjection) { matHelper.SetDoRadialProjection(fTargetZ); }
- matHelper.SetNraysPerDim(fMatBudgetNrays);
-
- double zLast = fTargetZ + 1.; // some gap (+1cm) to skip the target material
-
- for (unsigned int ist = 0; ist < vStations.size(); ist++) {
- auto& station = vStations[ist];
- double z1 = station.GetZdouble() + station.GetZthickness()[0] / 2.;
- double z2 = z1;
- if (ist < vStations.size() - 1) {
- // split materials between the stations at the middle
- auto& stationNext = vStations[ist + 1];
- z2 = stationNext.GetZdouble() - stationNext.GetZthickness()[0] / 2.;
- }
- double zNew = 0.5 * (z1 + z2);
-
- double maxXY = 1.3 * station.GetRmax()[0];
- //double maxXY = 80;
- // calculate n bins from the minimal pitch
- int nBins = static_cast<int>(std::ceil(2. * maxXY / fMatBudgetPitch));
- if (nBins < 1) { LOG(fatal) << " material nBins " << nBins << " is not positive, something is wrong"; }
- if (nBins > fMatBudgetNbins) { nBins = fMatBudgetNbins; }
-
- L1Material matBudget = matHelper.GenerateMaterialMap(station.GetZdouble(), zLast, zNew, maxXY, nBins);
-
- station.SetMaterialMap(matBudget);
-
- LOG(info) << "- Generated material map for tracking station " << ist << " at z = " << station.GetZdouble()
- << " cm."
- << " Material is collected between z = " << zLast << " and z = " << zNew;
-
- zLast = zNew;
- }
-
- auto timerEnd = std::chrono::high_resolution_clock::now();
- double materialGenerationTime = (double) (std::chrono::duration<double>(timerEnd - timerStart).count());
- LOG(info) << "- Generating material maps tooks " << materialGenerationTime << " seconds";
- }
-
- // ** Put stations to the init manager
-
- for (unsigned int ist = 0; ist < vStations.size(); ist++) {
- auto& station = vStations[ist];
- fInitManager.AddStation(station);
- }
// ****************************************
// ** **
@@ -968,14 +909,16 @@ InitStatus CbmL1::Init()
DumpMaterialToFile("L1material.root");
+ // Initialize counters
+ fEventNo = 0;
+ fNofRecoTracks = 0;
+
return kSUCCESS;
}
void CbmL1::Reconstruct(CbmEvent* event)
{
- static int nevent = 0;
-
fvSelectedMcEvents.clear();
int bestMcFile = -1;
@@ -1013,7 +956,7 @@ void CbmL1::Reconstruct(CbmEvent* event)
}
}
- if (fVerbose > 1) { cout << "\nCbmL1::Exec event " << ++nevent << " ...\n\n"; }
+ if (fVerbose > 1) { cout << "\nCbmL1::Exec event " << fEventNo << " ...\n\n"; }
#ifdef _OPENMP
omp_set_num_threads(1);
#endif
@@ -1021,6 +964,13 @@ void CbmL1::Reconstruct(CbmEvent* event)
// ----- Read data from branches and send data from IODataManager to L1Algo ----------------------------------------
ReadEvent(event);
+ //if (fPerformance) {
+ // LOG(info) << "**** DEBUG: BEGIN";
+ // LOG(info) << fvMCPoints[0].ToString(10, true); // header
+ // for (const auto& point: fvMCPoints) { LOG(info) << point.ToString(10); }
+ // LOG(info) << "**** DEBUG: END";
+ //}
+
// Material monitoring: mark active areas
{
@@ -1096,6 +1046,8 @@ void CbmL1::Reconstruct(CbmEvent* event)
fvRecoTracks.push_back(t);
}
+ fNofRecoTracks += static_cast<int>(fpAlgo->fRecoTracks.size());
+
LOG(debug) << "CA Track Finder: " << fpAlgo->fCaRecoTime << " s/sub-ts";
@@ -1131,6 +1083,8 @@ void CbmL1::Reconstruct(CbmEvent* event)
if ((symbol == 'e') || (symbol == 'q')) exit(0);
} while (symbol != '\n');
}
+
+ ++fEventNo;
}
// ----- Finish CbmStsFitPerformanceTask task -----------------------------
@@ -1138,15 +1092,23 @@ void CbmL1::Finish()
{
// monitor the material
- {
- LOG(info) << "\033[31;1m-------------------- L1 material -----------------------------\033[0m";
+ LOG(info) << "\033[31;1m ****************************\033[0m";
+ LOG(info) << "\033[31;1m ** CA Tracking monitore **\033[0m";
+ LOG(info) << "\033[31;1m ****************************\033[0m";
- for (int i = 0; i < fNStations; i++) {
- LOG(info) << fMaterialMonitor[i].ToString();
- }
- LOG(info) << "\033[31;1m--------------------------------------------------------------\033[0m";
+ LOG(info) << "\033[31;1m ----- Material budget -------------------------------------- \033[0m";
+ for (int i = 0; i < fNStations; i++) {
+ LOG(info) << fMaterialMonitor[i].ToString();
}
+ LOG(info) << "\033[31;1m -------------------------------------------------------------\033[0m";
+
+ // monitor of the reconstructed tracks
+ LOG(info) << "\033[31;1m ----- Counters ----------------------------------------------\033[0m";
+ LOG(info) << "\tNumber of analyzed events: " << fEventNo;
+ LOG(info) << "\tNumber of reconstructed tracks: " << fNofRecoTracks;
+ LOG(info) << "\tNumber of reconstructed tracks per event/TS: " << static_cast<double>(fNofRecoTracks) / fEventNo;
+ LOG(info) << "\033[31;1m -------------------------------------------------------------\033[0m";
TDirectory* curr = gDirectory;
TFile* currentFile = gFile;
@@ -1188,6 +1150,8 @@ void CbmL1::Finish()
}
+// ---------------------------------------------------------------------------------------------------------------------
+//
void CbmL1::writedir2current(TObject* obj)
{
if (!obj->IsFolder()) obj->Write();
@@ -1203,8 +1167,57 @@ void CbmL1::writedir2current(TObject* obj)
}
}
-/// ----- Ideal Tracking -----------------------------
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void CbmL1::GenerateMaterialMaps()
+{
+ LOG(info) << "Generating material maps...";
+ auto timerStart = std::chrono::high_resolution_clock::now();
+
+ ca::tools::MaterialHelper matHelper;
+ matHelper.SetSafeMaterialInitialization(fDoSafeMaterialInitialization);
+
+ if (!fMatBudgetParallelProjection) { matHelper.SetDoRadialProjection(fTargetZ); }
+ matHelper.SetNraysPerDim(fMatBudgetNrays);
+
+ double zLast = fTargetZ + 1.; // some gap (+1cm) to skip the target material
+
+ auto& vStations = fInitManager.GetStationInfo();
+ for (unsigned int ist = 0; ist < vStations.size(); ist++) {
+ auto& station = vStations[ist];
+ double z1 = station.GetZdouble() + station.GetZthickness()[0] / 2.;
+ double z2 = z1;
+ if (ist < vStations.size() - 1) {
+ // split materials between the stations at the middle
+ auto& stationNext = vStations[ist + 1];
+ z2 = stationNext.GetZdouble() - stationNext.GetZthickness()[0] / 2.;
+ }
+ double zNew = 0.5 * (z1 + z2);
+
+ double maxXY = 1.3 * station.GetRmax()[0];
+ //double maxXY = 80;
+ // calculate n bins from the minimal pitch
+ int nBins = static_cast<int>(std::ceil(2. * maxXY / fMatBudgetPitch));
+ if (nBins < 1) { LOG(fatal) << " material nBins " << nBins << " is not positive, something is wrong"; }
+ if (nBins > fMatBudgetNbins) { nBins = fMatBudgetNbins; }
+
+ L1Material matBudget = matHelper.GenerateMaterialMap(station.GetZdouble(), zLast, zNew, maxXY, nBins);
+
+ station.SetMaterialMap(matBudget);
+
+ LOG(info) << "Generated material map for tracking station " << ist << " at z = " << station.GetZdouble() << " cm."
+ << " Material is collected between z = " << zLast << " and z = " << zNew;
+ zLast = zNew;
+ }
+
+ auto timerEnd = std::chrono::high_resolution_clock::now();
+ double materialGenerationTime = (double) (std::chrono::duration<double>(timerEnd - timerStart).count());
+ LOG(info) << "Generating material maps... done! (it took " << materialGenerationTime << " s)";
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
void CbmL1::IdealTrackFinder()
{
fpAlgo->fRecoTracks.clear();
diff --git a/reco/L1/CbmL1.h b/reco/L1/CbmL1.h
index df9317cf0f5d471a797408d96d0422669662e811..e2c580c1a732982ab45b3a3069fdd090e6aea527 100644
--- a/reco/L1/CbmL1.h
+++ b/reco/L1/CbmL1.h
@@ -336,13 +336,17 @@ private:
float xMin, xMax;
};
- /// Runs ideal track finder: copies all MC-tracks into reconstructed tracks
- void IdealTrackFinder();
+ /// @brief Generates material maps
+ void GenerateMaterialMaps();
- /// Read information about hits, mcPoints and mcTracks into L1 classes
+ /// @brief Runs ideal track finder: copies all MC-tracks into reconstructed tracks
+ void IdealTrackFinder();
- /// Repacks data from the external TClonesArray objects to the internal L1 arrays
- /// \param event Pointer to the current CbmEvent object
+ /// @brief Reads input data in the event
+ /// @param event Pointer to the current CbmEvent object
+ ///
+ /// Reads information about hits, mcPoints and mcTracks into L1 classes. Repacks data from the external data branches
+ /// to the internal data arrays.
void ReadEvent(CbmEvent* event = NULL);
/// Converts data from generic FairMCPoint based class to the CbmL1MCPoint
@@ -636,6 +640,9 @@ private:
L1Vector<CbmL1MCTrack> fvMCTracks = {"CbmL1::fvMCTracks"}; ///< Array of MC tracks
L1Vector<int> fvHitPointIndexes = {"CbmL1::fvHitPointIndexes"}; ///< Indexes of MC points vs. hit index
+ int fEventNo = 0; ///< Current number of event/TS
+ int fNofRecoTracks = 0; ///< Total number of reconstructed tracks
+
public:
// ** Repacked input data **
diff --git a/reco/L1/CbmL1DetectorID.h b/reco/L1/CbmL1DetectorID.h
index 90c278fc4a1f6971d965801373b8743b12b94c96..2ae838bd02d925627b591b8dcd75682635a39bd1 100644
--- a/reco/L1/CbmL1DetectorID.h
+++ b/reco/L1/CbmL1DetectorID.h
@@ -19,9 +19,10 @@
/// 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
+/// NOTE: The enumeration must not contain jumps in ordering and the first entry must be equal 0
enum class L1DetectorID
{
- kMvd,
+ kMvd = 0,
kSts,
kMuch,
kTrd,
diff --git a/reco/L1/CbmL1Performance.cxx b/reco/L1/CbmL1Performance.cxx
index 9f02e032d784f6bfdf2f1c040958041079b96ac4..6f44303449d1585daaa5e7df284d72ad5c042c2f 100644
--- a/reco/L1/CbmL1Performance.cxx
+++ b/reco/L1/CbmL1Performance.cxx
@@ -75,7 +75,6 @@ void CbmL1::TrackMatch()
for (int iTrk = 0; iTrk < int(fvMCTracks.size()); ++iTrk) {
assert(iTrk == fvMCTracks[iTrk].ID);
}
-
// ***** DEBUG: End
// fill pMCTrackMap
diff --git a/reco/L1/CbmL1ReadEvent.cxx b/reco/L1/CbmL1ReadEvent.cxx
index 42810753fad53423f84c9e7cffd232622d525d7f..1453dc42cd170c110f6098e807d4017baa51e336 100644
--- a/reco/L1/CbmL1ReadEvent.cxx
+++ b/reco/L1/CbmL1ReadEvent.cxx
@@ -26,6 +26,7 @@
#include "CbmMuchGeoScheme.h"
#include "CbmMuchPixelHit.h"
#include "CbmMuchPoint.h"
+#include "CbmMuchTrackingInterface.h"
#include "CbmStsAddress.h"
#include "CbmStsCluster.h"
#include "CbmStsDigi.h"
@@ -319,12 +320,9 @@ void CbmL1::ReadEvent(CbmEvent* event)
for (Int_t iMC = 0; iMC < nMC; iMC++) {
CbmL1MCPoint MC;
if (!ReadMCPoint(&MC, iMC, iFile, iEvent, 2)) {
- MC.iStation = -1;
- const L1Station* sta = fpAlgo->GetParameters()->GetStations().begin();
- for (Int_t iSt = 0; iSt < fNMuchStationsGeom; iSt++) {
- int iStActive = fpAlgo->GetParameters()->GetStationIndexActive(iSt, L1DetectorID::kMuch);
- if (MC.z > sta[iStActive].fZ[0] - 2.5) { MC.iStation = iStActive; }
- }
+ auto* pPoint = static_cast<CbmMuchPoint*>(fpMuchPoints->Get(iFile, iEvent, iMC));
+ int iStLoc = CbmMuchTrackingInterface::Instance()->GetTrackingStationIndex(pPoint->GetDetectorID());
+ MC.iStation = fpAlgo->GetParameters()->GetStationIndexActive(iStLoc, L1DetectorID::kMuch);
if (MC.iStation < 0) { continue; } // Reject MC points from inactive stations;
auto itTrack = fmMCTracksLinksMap.find(CbmL1LinkKey(MC.ID, iEvent, iFile));
assert(itTrack != fmMCTracksLinksMap.cend());
diff --git a/reco/L1/L1Algo/L1BaseStationInfo.cxx b/reco/L1/L1Algo/L1BaseStationInfo.cxx
index ae211888c0cf3956659e90f020b2cd728acdade0..667badd8e45c9bc2ebc8ce188cd01abbca751d3f 100644
--- a/reco/L1/L1Algo/L1BaseStationInfo.cxx
+++ b/reco/L1/L1Algo/L1BaseStationInfo.cxx
@@ -138,6 +138,7 @@ void L1BaseStationInfo::Reset()
//
const L1Station& L1BaseStationInfo::GetL1Station() const
{
+ // TODO: replace L1MASSERT with throw std::logic_error
std::stringstream aStream;
aStream << "L1BaseStationInfo::GetL1Station: attempt to get an L1Staion object from uninitialized L1BaseStation with "
<< "stationID = " << fStationID << " and detectorID = " << static_cast<int>(fDetectorID);
diff --git a/reco/L1/L1Algo/L1BaseStationInfo.h b/reco/L1/L1Algo/L1BaseStationInfo.h
index 316c590d823011bc84cc3af03d10ae4935a05478..4045f81ac5b0037767b9a7bef5abae37ca06079d 100644
--- a/reco/L1/L1Algo/L1BaseStationInfo.h
+++ b/reco/L1/L1Algo/L1BaseStationInfo.h
@@ -91,10 +91,6 @@ public:
bool operator<(const L1BaseStationInfo& right) const
{
return (GetZdouble() < right.GetZdouble());
-
- /// This way sorting is carried out first by fDetectorID,
- /// and then by fStationID. The detectorID must be defined in L1DetectorID enum class
- // return fDetectorID != right.fDetectorID ? fDetectorID < right.fDetectorID : fStationID < right.fStationID;
}
/// Gets detector ID
@@ -204,7 +200,7 @@ public:
void SetRmin(double inRmin);
/// Sets station ID
- void SetStationID(int inID);
+ [[deprecated("Please, use constructor to set station ID")]] void SetStationID(int inID);
/// Sets type of station
void SetStationType(int inType); // TODO: this is a temporary solution (S.Zh.)
@@ -239,7 +235,7 @@ public:
private:
L1DetectorID fDetectorID {static_cast<L1DetectorID>(0)}; ///< Detector ID
- int fStationID {-1}; ///< Station ID
+ int fStationID {-1}; ///< Local ID of a station
bool fTrackingStatus {false}; ///< Tracking status: true - station is used for tracking
double fXmax {0}; ///< Maximum distance between station center and its edge in x direction
double fYmax {0}; ///< Maximum distance between station center and its edge in y direction
diff --git a/reco/L1/L1Algo/L1Def.h b/reco/L1/L1Algo/L1Def.h
index b372c25cabc22999546ec08512024f21b8c6735e..0e7da09ac99f4d5d67c4d21f645a7202572ba59e 100644
--- a/reco/L1/L1Algo/L1Def.h
+++ b/reco/L1/L1Algo/L1Def.h
@@ -54,8 +54,14 @@
std::cout << __FILE__ << ":" << __LINE__ << ": \033[01;38;5;208m" << (#expr) << "\033[0m = " << (expr) << "\n"
#define L1_SHOWF(msg) \
std::cout << "(!) " << __FILE__ << ":" << __LINE__ << ": \033[01;38;5;208m" << (#msg) << "\033[0m\n"
+#define L1_SHOWCONTAINER(cont) \
+ std::cout << __FILE__ << ":" << __LINE__ << ": \033[01;38;5;208m" << (#cont) << "\033[0m: "; \
+ std::for_each(cont.cbegin(), cont.cend(), [](const auto& el) { std::cout << el << ' '; }); \
+ std::cout << '\n';
+
#endif
+
// Prints function call
#if defined(__GNUC__)
#define L1_SHOWFN std::cout << "\033[1;32mCALL \033[1;33m" << __PRETTY_FUNCTION__ << "\033[0m\n"
diff --git a/reco/L1/L1Algo/L1IODataManager.cxx b/reco/L1/L1Algo/L1IODataManager.cxx
index 37913828d6ff099049c8188bfbac8c61c09d75cc..66d695fc590552099f21f7d9c13f4b08dff2b348 100644
--- a/reco/L1/L1Algo/L1IODataManager.cxx
+++ b/reco/L1/L1Algo/L1IODataManager.cxx
@@ -93,6 +93,7 @@ void L1IODataManager::InitData()
//std::cout << "N data streams: " << fInputData.fStreamStartIndices.size() << std::endl;
+ // TODO: SZh 14.08.2023: Move the max allowed number of streams to the L1Constants.h
if (fInputData.fStreamStartIndices.size() > 3000) {
LOG(warning) << "L1: unexpected order of input data: too many data streams!!! ";
fInputData.fStreamStartIndices.reduce(3000);
diff --git a/reco/L1/L1Algo/L1InitManager.cxx b/reco/L1/L1Algo/L1InitManager.cxx
index 745fc35d86394ff41ba32b904f63eba8bf9769a1..2bd6fcda99233bdb4a3c8de1a834e974d93324ae 100644
--- a/reco/L1/L1Algo/L1InitManager.cxx
+++ b/reco/L1/L1Algo/L1InitManager.cxx
@@ -27,61 +27,14 @@ using L1Constants::clrs::kRDb; // bold red log
//
void L1InitManager::AddStation(const L1BaseStationInfo& inStation)
{
- // Check if other fields were defined already
- // Active detector IDs
+ // TODO: SZh 15.08.2023: Replace L1MASSERT with throw logic_error
L1MASSERT(0, fInitController.GetFlag(EInitKey::kActiveDetectorIDs),
"Attempt to add a station info before the active detectors set had been initialized");
- // Number of stations check
- L1MASSERT(0, fInitController.GetFlag(EInitKey::kStationsNumberCrosscheck),
- "Attempt to add a station info before the numbers of stations for each detector had been initialized");
-
- // Field function
- L1MASSERT(0, fInitController.GetFlag(EInitKey::kFieldFunction),
- "Attempt to add a station info before the magnetic field function had been initialized");
-
- // Check activeness of this station type
- bool isStationActive =
- inStation.GetTrackingStatus() && fActiveDetectorIDs.find(inStation.GetDetectorID()) != fActiveDetectorIDs.end();
- if (isStationActive) {
- // initialize magnetic field slice
- L1BaseStationInfo inStationCopy = L1BaseStationInfo(inStation); // make a copy of station so it can be initialized
- inStationCopy.SetFieldFunction(fFieldFunction);
-
- // check, if material map is set
- if (!inStationCopy.GetInitController().GetFlag(L1BaseStationInfo::EInitKey::kThicknessMap)) {
- LOG(fatal) << "Station material map was not set for detectorID = "
- << static_cast<int>(inStationCopy.GetDetectorID()) << ", stationID = " << inStationCopy.GetStationID();
- }
-
- // Check station init
- LOG(debug) << "L1InitManager::AddStation:(original) L1BaseStationInfo " << inStation.GetInitController().ToString();
- LOG(debug) << "L1InitManager::AddStation:(copy) L1BaseStationInfo " << inStation.GetInitController().ToString();
- if (!inStationCopy.GetInitController().IsFinalized()) {
- LOG(fatal) << "Attempt to add an incompletely initialized station info object (detectorID = "
- << static_cast<int>(inStationCopy.GetDetectorID()) << ", stationID = " << inStationCopy.GetStationID()
- << ")";
- }
- // insert the station in a set
- auto insertionResult = fStationsInfo.insert(std::move(inStationCopy));
- if (!insertionResult.second) {
- LOG(fatal) << "Attempt to add a duplicated station info object (detectorID = "
- << static_cast<int>(inStationCopy.GetDetectorID()) << ", stationID = " << inStationCopy.GetStationID()
- << ")";
- }
-
- int index = fStationsInfo.size() - 1 + fParameters.fNstationsGeometryTotal - fParameters.fNstationsActiveTotal;
- fParameters.fActiveStationGlobalIDs[index] = fStationsInfo.size() - 1;
- }
- else {
- int index = fStationsInfo.size() + fParameters.fNstationsGeometryTotal - fParameters.fNstationsActiveTotal;
- fParameters.fActiveStationGlobalIDs[index] = -1;
- fParameters.fNstationsActive[static_cast<L1DetectorID_t>(inStation.GetDetectorID())]--;
- fParameters.fNstationsActiveTotal--;
+ // Check, if the detector subsystem for this station is active
+ if (fActiveDetectorIDs.find(inStation.GetDetectorID()) != fActiveDetectorIDs.end()) {
+ fvStationInfo.push_back(inStation);
}
- LOG(debug) << "L1InitManager: adding a station with stationID = " << inStation.GetStationID()
- << " and detectorID = " << static_cast<int>(inStation.GetDetectorID())
- << ". Is active: " << isStationActive;
}
// ----------------------------------------------------------------------------------------------------------------------
@@ -97,17 +50,17 @@ void L1InitManager::CheckInit()
void L1InitManager::ClearSetupInfo()
{
// Clear stations set and a thickness map
- fStationsInfo.clear();
+ fvStationInfo.clear();
std::fill(fParameters.fThickMap.begin(), fParameters.fThickMap.end(), L1Material());
fInitController.SetFlag(EInitKey::kStationsInfo, false);
// Set number of stations do default values
- std::fill(fParameters.fNstationsGeometry.begin(), fParameters.fNstationsGeometry.end(), 0);
- std::fill(fParameters.fNstationsActive.begin(), fParameters.fNstationsActive.end(), 0);
- std::fill(fParameters.fActiveStationGlobalIDs.begin(), fParameters.fActiveStationGlobalIDs.end(), 0);
- fParameters.fNstationsGeometryTotal = -1;
+ std::fill(fParameters.fvFirstGeoId.begin(), fParameters.fvFirstGeoId.end(), 0);
+ std::fill(fParameters.fvLocalToGeoIdMap.begin(), fParameters.fvLocalToGeoIdMap.end(), 0);
+ std::fill(fParameters.fvGeoToActiveMap.begin(), fParameters.fvGeoToActiveMap.end(), 0);
+ std::fill(fParameters.fvActiveToGeoMap.begin(), fParameters.fvActiveToGeoMap.end(), 0);
fParameters.fNstationsActiveTotal = -1;
- fInitController.SetFlag(EInitKey::kStationsNumberCrosscheck, false);
+ fInitController.SetFlag(EInitKey::kStationLayoutInitialized, false);
// Clear active detectors
fActiveDetectorIDs.clear();
@@ -174,6 +127,9 @@ bool L1InitManager::FormParametersContainer()
if (!fParameters.fDevIsParSearchWUsed) { fInitController.SetFlag(EInitKey::kSearchWindows, true); }
+ // Apply magnetic field to the station info objects
+ std::for_each(fvStationInfo.begin(), fvStationInfo.end(), [&](auto& st) { st.SetFieldFunction(fFieldFunction); });
+
// Check initialization
this->CheckInit();
@@ -183,20 +139,22 @@ bool L1InitManager::FormParametersContainer()
return false;
}
- // Form array of stations
- auto destinationArrayIterator = fParameters.fStations.begin();
- for (const auto& item : fStationsInfo) {
- *destinationArrayIterator = item.GetL1Station();
- ++destinationArrayIterator;
+ { // Form array of stations
+ auto destIt = fParameters.fStations.begin();
+ for (const auto& station : fvStationInfo) {
+ if (!station.GetTrackingStatus()) { continue; }
+ *destIt = station.GetL1Station();
+ ++destIt;
+ }
}
- // Form array of material map
- auto thickMapIt = fParameters.fThickMap.begin();
- for (auto it = fStationsInfo.begin(); it != fStationsInfo.end(); ++it) {
- auto node = fStationsInfo.extract(it);
- *thickMapIt = std::move(node.value().TakeMaterialMap());
- fStationsInfo.insert(std::move(node));
- ++thickMapIt;
+ { // Form array of material map
+ auto destIt = fParameters.fThickMap.begin();
+ for (auto& station : fvStationInfo) {
+ if (!station.GetTrackingStatus()) { continue; }
+ *destIt = std::move(station.TakeMaterialMap());
+ ++destIt;
+ }
}
// Check the consistency of the parameters object. If object inconsistent, it throws std::logic_error
@@ -211,6 +169,81 @@ bool L1InitManager::FormParametersContainer()
return true;
}
+// ----------------------------------------------------------------------------------------------------------------------
+//
+int L1InitManager::GetNstationsActive() const
+{
+ LOG_IF(fatal, !fInitController.GetFlag(EInitKey::kStationLayoutInitialized))
+ << "L1InitManager: number of active stations cannot be accessed until the station layout is initialized";
+ return fParameters.GetNstationsActive();
+}
+
+// ----------------------------------------------------------------------------------------------------------------------
+//
+int L1InitManager::GetNstationsActive(L1DetectorID detectorID) const
+{
+ LOG_IF(fatal, !fInitController.GetFlag(EInitKey::kStationLayoutInitialized))
+ << "L1InitManager: number of active stations cannot be accessed until the station layout is initialized";
+ return fParameters.GetNstationsActive(detectorID);
+}
+
+// ----------------------------------------------------------------------------------------------------------------------
+//
+int L1InitManager::GetNstationsGeometry() const
+{
+ LOG_IF(fatal, !fInitController.GetFlag(EInitKey::kStationLayoutInitialized))
+ << "L1InitManager: number of geometry stations cannot be accessed until the station layout is initialized";
+ return fParameters.GetNstationsGeometry();
+}
+
+// ----------------------------------------------------------------------------------------------------------------------
+//
+int L1InitManager::GetNstationsGeometry(L1DetectorID detectorID) const
+{
+ LOG_IF(fatal, !fInitController.GetFlag(EInitKey::kStationLayoutInitialized))
+ << "L1InitManager: number of geometry stations cannot be accessed until the station layout is initialized";
+ return fParameters.GetNstationsGeometry(detectorID);
+}
+
+// ----------------------------------------------------------------------------------------------------------------------
+//
+std::vector<L1BaseStationInfo>& L1InitManager::GetStationInfo()
+{
+ LOG_IF(fatal, !fInitController.GetFlag(EInitKey::kStationLayoutInitialized))
+ << "L1InitManager: station info container cannot be accessed until the station layout is initialized";
+ return fvStationInfo;
+}
+
+
+// ----------------------------------------------------------------------------------------------------------------------
+//
+void L1InitManager::InitStationLayout()
+{
+ std::sort(fvStationInfo.begin(), fvStationInfo.end());
+
+ for (const auto& aStation : fvStationInfo) {
+ ++fParameters.fvFirstGeoId[static_cast<int>(aStation.GetDetectorID()) + 1];
+ }
+ for (int iDet = 1; iDet < static_cast<int>(fParameters.fvFirstGeoId.size()); ++iDet) {
+ fParameters.fvFirstGeoId[iDet + 1] += fParameters.fvFirstGeoId[iDet];
+ }
+
+ fParameters.fNstationsActiveTotal = 0;
+ for (int iStGeo = 0; iStGeo < static_cast<int>(fvStationInfo.size()); ++iStGeo) {
+ const auto& aStation = fvStationInfo[iStGeo];
+ int iDet = static_cast<int>(aStation.GetDetectorID());
+ int iStLocal = aStation.GetStationID();
+ // Fill geo -> local map
+ fParameters.fvLocalToGeoIdMap[fParameters.fvFirstGeoId[iDet] + iStLocal] = iStGeo;
+ // Fill geo <-> active map
+ int iStActive = aStation.GetTrackingStatus() ? fParameters.fNstationsActiveTotal++ : -1;
+ fParameters.fvGeoToActiveMap[iStGeo] = iStActive;
+ if (iStActive > -1) { fParameters.fvActiveToGeoMap[iStActive] = iStGeo; }
+ }
+
+ fInitController.SetFlag(EInitKey::kStationLayoutInitialized, true);
+}
+
// ----------------------------------------------------------------------------------------------------------------------
//
void L1InitManager::InitTargetField(double zStep)
@@ -403,47 +436,6 @@ void L1InitManager::SetMomentumCutOff(float momentumCutOff)
fInitController.SetFlag(EInitKey::kMomentumCutOff);
}
-// ----------------------------------------------------------------------------------------------------------------------
-//
-void L1InitManager::SetNstations(L1DetectorID detectorID, int nStations)
-{
- L1MASSERT(0, fInitController.GetFlag(EInitKey::kActiveDetectorIDs),
- "Attempt to set crosscheck number of stations before the active detetors set had been initialized");
-
- // NOTE: We add and check only those detectors which will be active (?)
- // For INACTIVE detectors the initialization code for it inside CbmL1/BmnL1 can (and must) be still in,
- // but it will be ignored inside L1InitManager.
- if (fActiveDetectorIDs.find(detectorID) != fActiveDetectorIDs.end()) {
- if (nStations) {
- fParameters.fNstationsGeometry[static_cast<L1DetectorID_t>(detectorID)] = nStations;
- fParameters.fNstationsActive[static_cast<L1DetectorID_t>(detectorID)] = nStations;
- }
- else {
- // TODO: Probably it is better to replace fatal with warn and remove the detectorID from active detectors (S.Zharko)
- LOG(fatal) << "L1InitManager::SetNstations: attempt to initialize zero stations for active detector: "
- << static_cast<L1DetectorID_t>(detectorID);
- }
- }
-
- // Check if all the station numbers for active detectors are initialized now:
- LOG(debug) << "L1InitManager::SetNstations called for detectorID = " << static_cast<int>(detectorID);
- if (!fInitController.GetFlag(EInitKey::kStationsNumberCrosscheck)) {
- bool ifInitialized = true;
- for (auto item : fActiveDetectorIDs) {
- if (fParameters.fNstationsGeometry[static_cast<L1DetectorID_t>(item)] == 0) {
- ifInitialized = false;
- break;
- }
- }
- fInitController.SetFlag(EInitKey::kStationsNumberCrosscheck, ifInitialized);
- }
- if (fInitController.GetFlag(EInitKey::kStationsNumberCrosscheck)) {
- fParameters.fNstationsGeometryTotal =
- std::accumulate(fParameters.fNstationsGeometry.begin(), fParameters.fNstationsGeometry.end(), 0);
- fParameters.fNstationsActiveTotal = fParameters.fNstationsGeometryTotal;
- }
-}
-
// ----------------------------------------------------------------------------------------------------------------------
//
void L1InitManager::SetTargetPosition(double x, double y, double z)
@@ -522,31 +514,14 @@ void L1InitManager::CheckStationsInfoInit()
{
bool ifInitPassed = true;
if (!fInitController.GetFlag(EInitKey::kStationsInfo)) {
- //
- // 1) Check numbers of stations passed
- //
- // loop over active detectors
- for (auto itemDetector : fActiveDetectorIDs) {
- auto selectDetector = [&itemDetector](const L1BaseStationInfo& station) {
- return station.GetDetectorID() == itemDetector;
- };
- int nStationsExpected = GetNstationsActive(itemDetector);
- int nStations = std::count_if(fStationsInfo.begin(), fStationsInfo.end(), selectDetector);
- if (nStations != nStationsExpected) {
- LOG(error) << "L1InitManager::CheckStationsInfoInit: Incorrect number of L1BaseStationInfo objects passed"
- << " to the L1Manager for L1DetectorID = " << static_cast<int>(itemDetector) << ": " << nStations
- << " of " << nStationsExpected << " expected";
- ifInitPassed = false;
- }
- } // loop over active detectors: end
- L1MASSERT(0, ifInitPassed, "Station info initialization failed");
-
- //
- // 2) Check for maximum allowed number of stations
- //
- int nStationsTotal = fParameters.fNstationsGeometryTotal;
- if (nStationsTotal > L1Constants::size::kMaxNstations) {
- LOG(fatal) << "Actual total number of registered stations in geometry (" << nStationsTotal
+ // (1) Check the stations themselves
+ bool bStationsFinalized = std::all_of(fvStationInfo.begin(), fvStationInfo.end(),
+ [](const auto& st) { return st.GetInitController().IsFinalized(); });
+ if (!bStationsFinalized) { LOG(fatal) << "At least one of the L1BaseStationInfo objects is not finalized"; }
+
+ // (2) Check for maximum allowed number of stations
+ if (fParameters.GetNstationsGeometry() > L1Constants::size::kMaxNstations) {
+ LOG(fatal) << "Actual total number of registered stations in geometry (" << fParameters.GetNstationsGeometry()
<< ") is larger then possible (" << L1Constants::size::kMaxNstations
<< "). Please, select another set of active tracking detectors or recompile the code with enlarged"
<< " L1Constants::size::kMaxNstations value";
diff --git a/reco/L1/L1Algo/L1InitManager.h b/reco/L1/L1Algo/L1InitManager.h
index e85cc484ab3d252f128e7c6f6e659ec9c65ccafe..4629b5090af0ef52fdea1afa7ec5a713c4362ddf 100644
--- a/reco/L1/L1Algo/L1InitManager.h
+++ b/reco/L1/L1Algo/L1InitManager.h
@@ -50,6 +50,7 @@ using L1DetectorID_t = std::underlying_type_t<L1DetectorID>;
/// kTof,
/// kTrd
/// };
+/// For CBM see implementation in CbmL1DetectorID.h
///
/// 1. Get a pointer to the :L1InitManager field of L1Algo:
///
@@ -75,17 +76,17 @@ private:
// NOTE: Please, keep the numbers of the enumeration items in the existing order: it helps to debug the
// initialization with this->GetObjectInitController().ToString() method call (S.Zharko)
kActiveDetectorIDs, ///< 0) If the detector sequence is set
- kStationsNumberCrosscheck, ///< 1) If the crosscheck station numbers were setup
- kFieldFunction, ///< 2) If magnetic field getter function is set
- kTargetPos, ///< 3) If target position was defined
- kPrimaryVertexField, ///< 4) If magnetic field value and region defined at primary vertex
- kStationsInfo, ///< 5) If all the planned stations were added to the manager
- kCAIterationsNumberCrosscheck, ///< 6) If the number of CA track finder is initialized
- kCAIterations, ///< 7) If the CA track finder iterations were initialized
- kSearchWindows, ///< 8) If the hit search windows were initialized
- kTrackingLevel, ///< 9)
- kGhostSuppression, ///< 10)
- kMomentumCutOff, ///< 11)
+ kFieldFunction, ///< 1) If magnetic field getter function is set
+ kTargetPos, ///< 2) If target position was defined
+ kPrimaryVertexField, ///< 3) If magnetic field value and region defined at primary vertex
+ kStationsInfo, ///< 4) If all the planned stations were added to the manager
+ kCAIterationsNumberCrosscheck, ///< 5) If the number of CA track finder is initialized
+ kCAIterations, ///< 6) If the CA track finder iterations were initialized
+ kSearchWindows, ///< 7) If the hit search windows were initialized
+ kTrackingLevel, ///< 8)
+ kGhostSuppression, ///< 9)
+ kMomentumCutOff, ///< 10)
+ kStationLayoutInitialized, ///< 11) If stations layout is initialized
kEnd ///< 12) [technical] number of entries in the enumeration
};
@@ -114,15 +115,11 @@ public:
L1InitManager& operator=(L1InitManager&& /*other*/) = delete;
- /// Adds another station of a given type using reference to a L1BaseStationInfo object
+ /// @brief Adds a tracking station to the geometry
+ ///
+ /// @note The added station stays uninitialized until the L1Parameters object is formed
void AddStation(const L1BaseStationInfo& station);
- /// Adds another station of a given type using pointer to a L1BaseStationInfo object
- void AddStation(const L1BaseStationInfo* pStation) { AddStation(*pStation); }
-
- /// Adds another station of a given type using std::unique_ptr-wraped pointer to L1BaseStationInfo
- void AddStation(const std::unique_ptr<L1BaseStationInfo>& puStation) { AddStation(*puStation); }
-
/// Provides final checks of large fields initialization calling Check"Object"Init() private methods,
/// must be called in the beginning of L1Algo::Init()
void CheckInit();
@@ -152,30 +149,34 @@ public:
const InitController_t& GetInitController() const { return fInitController; }
/// Gets total number of active stations
- int GetNstationsActive() const { return fParameters.fNstationsActiveTotal; }
+ [[deprecated]] int GetNstationsActive() const;
/// Gets number of active stations for given detector ID
- int GetNstationsActive(L1DetectorID detectorID) const
- {
- return fParameters.fNstationsActive[static_cast<L1DetectorID_t>(detectorID)];
- }
+ [[deprecated]] int GetNstationsActive(L1DetectorID detectorID) const;
/// Gets total number of stations, provided by setup geometry
- int GetNstationsGeometry() const { return fParameters.fNstationsGeometryTotal; }
+ [[deprecated]] int GetNstationsGeometry() const;
/// Gets number of stations, provided by setup geometry for given detector ID
- int GetNstationsGeometry(L1DetectorID detectorID) const
- {
- return fParameters.fNstationsGeometry[static_cast<L1DetectorID_t>(detectorID)];
- }
+ [[deprecated]] int GetNstationsGeometry(L1DetectorID detectorID) const;
/// Gets a name of the output configuration file
const std::string& GetOutputConfigName() const { return fConfigOutputName; }
+ /// @brief Gets a reference to the stations array
+ std::vector<L1BaseStationInfo>& GetStationInfo();
+
/// Gets tracking level
int GetTrackingLevel() const { return fParameters.fTrackingLevel; }
- /// Calculates L1FieldValue and L1FieldReference values for a selected step in z-axis from the target position
+ /// @brief Initializes station layout
+ ///
+ /// This function is to be called after all the tracking stations (L1BaseStationInfo objects) are added to the
+ /// L1InitManager instance. After the initialization a vector of the tracking stations sorted by z-positions is
+ /// available for modifications.
+ void InitStationLayout();
+
+ /// @brief Calculates L1FieldValue and L1FieldReference values for a selected step in z-axis from the target position
/// \param zStep step between nodal points
// TODO: Consider possibility for linear approximation (S.Zh.)
void InitTargetField(double zStep);
@@ -217,9 +218,6 @@ public:
///
void SetMomentumCutOff(float momentumCutOff);
- /// Sets a number of actual stations for a particular tracking detector ID to provide initialization cross-check
- void SetNstations(L1DetectorID detectorID, int nStations);
-
/// Sets a name of the output configuration file. The output file is created from the fields, saved in the resulted
/// L1Parameters object
/// \param filename Name of the output L1 parameters configuration
@@ -315,7 +313,7 @@ private:
// * Stations related fields *
- std::set<L1BaseStationInfo> fStationsInfo {}; ///< Set of L1BaseStationInfo objects
+ std::vector<L1BaseStationInfo> fvStationInfo {}; ///< Vector of L1BaseStationInfo objects (active + inactive)
/// A function which returns magnetic field vector B in a radius-vector xyz
L1FieldFunction_t fFieldFunction {[](const double (&)[3], double (&)[3]) {}};
diff --git a/reco/L1/L1Algo/L1Parameters.cxx b/reco/L1/L1Algo/L1Parameters.cxx
index 396c09b151a59ec9466bb4defa63510a353b853b..9ad3136445dda2681dc08c90da62d0d9b0177275 100644
--- a/reco/L1/L1Algo/L1Parameters.cxx
+++ b/reco/L1/L1Algo/L1Parameters.cxx
@@ -13,37 +13,18 @@
#include <iomanip>
-//----------------------------------------------------------------------------------------------------------------------
+// ---------------------------------------------------------------------------------------------------------------------
//
L1Parameters::L1Parameters()
{
- fActiveStationGlobalIDs.fill(-1); // by default, all stations are inactive, thus all the IDs must be -1
+ fvGeoToActiveMap.fill(-1); // by default, all stations are inactive, thus all the IDs must be -1
}
-//----------------------------------------------------------------------------------------------------------------------
+// ---------------------------------------------------------------------------------------------------------------------
//
void L1Parameters::CheckConsistency() const
{
LOG(info) << "Consistency test for L1 parameters object... ";
- /*
- * Arrays containing numbers of stations
- *
- * In the fNstationsActive and fNstationsGeometry array first L1Constants::size::kMaxNdetectors elements are the numbers of stations
- * in particular detector subsystem. The last element in both arrays corresponds to the total number of stations (geometry or
- * active). This fact applies restriction on the arrays: the last element must be equal to the sum of all previous elements.
- *
- */
-
- if (std::accumulate(fNstationsGeometry.cbegin(), fNstationsGeometry.cend(), 0) != fNstationsGeometryTotal) {
- throw std::logic_error("L1Parameters: invalid object condition: total number of stations provided by geometry "
- "differs from the sum of the station numbers for individual detector subsystems");
- };
-
- if (std::accumulate(fNstationsActive.cbegin(), fNstationsActive.cend(), 0) != fNstationsActiveTotal) {
- throw std::logic_error("L1Parameters: invalid object condition: total number of stations active in tracking "
- "differs from the sum of the station numbers for individual detector subsystems");
- };
-
/*
* Array of active station IDs
*
@@ -53,8 +34,7 @@ void L1Parameters::CheckConsistency() const
*/
auto filterInactiveStationIDs = [](int x) { return x != -1; };
- int nStationsCheck =
- std::count_if(fActiveStationGlobalIDs.cbegin(), fActiveStationGlobalIDs.cend(), filterInactiveStationIDs);
+ int nStationsCheck = std::count_if(fvGeoToActiveMap.cbegin(), fvGeoToActiveMap.cend(), filterInactiveStationIDs);
if (nStationsCheck != fNstationsActiveTotal) {
std::stringstream msg;
msg << "L1Parameters: invalid object condition: array of active station IDs is not consistent "
@@ -68,8 +48,7 @@ void L1Parameters::CheckConsistency() const
std::vector<int> idsCheck(
nStationsCheck); // temporary vector containing a sequence of active station IDs without -1 elements
- std::copy_if(fActiveStationGlobalIDs.cbegin(), fActiveStationGlobalIDs.cend(), idsCheck.begin(),
- filterInactiveStationIDs);
+ std::copy_if(fvGeoToActiveMap.cbegin(), fvGeoToActiveMap.cend(), idsCheck.begin(), filterInactiveStationIDs);
bool isStationIDsOk = true;
for (int id = 0; id < nStationsCheck; ++id) {
isStationIDsOk = isStationIDsOk && idsCheck[id] == id;
@@ -78,7 +57,7 @@ void L1Parameters::CheckConsistency() const
std::stringstream msg;
msg << "L1Parameters: invalid object condition: array of active station IDs is not a gapless subset "
<< "of integer numbers starting from 0:\n\t";
- for (auto id : fActiveStationGlobalIDs) {
+ for (auto id : fvGeoToActiveMap) {
msg << std::setw(3) << std::setfill(' ') << id << ' ';
}
throw std::logic_error(msg.str());
@@ -207,11 +186,23 @@ void L1Parameters::CheckConsistency() const
LOG(info) << "Consistency test for L1 parameters object... \033[1;32mpassed\033[0m";
}
-//----------------------------------------------------------------------------------------------------------------------
+// ---------------------------------------------------------------------------------------------------------------------
+//
+int L1Parameters::GetNstationsActive(L1DetectorID detectorID) const
+{
+ int nStations = 0;
+ for (int iStLoc = 0; iStLoc < this->GetNstationsGeometry(detectorID); ++iStLoc) {
+ int iStActive = this->GetStationIndexActive(iStLoc, detectorID);
+ if (iStActive > -1) { ++nStations; }
+ }
+ return nStations;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
//
void L1Parameters::Print(int /*verbosityLevel*/) const { LOG(info) << ToString(); }
-//----------------------------------------------------------------------------------------------------------------------
+// ---------------------------------------------------------------------------------------------------------------------
//
std::string L1Parameters::ToString(int verbosity, int indentLevel) const
{
@@ -252,20 +243,20 @@ std::string L1Parameters::ToString(int verbosity, int indentLevel) const
}
aStream << indent << indentChar << kCLb << "NUMBER OF STATIONS:\n" << kCL;
aStream << indent << indentChar << indentChar << "Number of stations (Geometry): ";
- for (int idx = 0; idx < int(fNstationsGeometry.size()); ++idx) {
- aStream << std::setw(2) << std::setfill(' ') << fNstationsGeometry[idx] << ' ';
+ for (int iDet = 0; iDet < L1Constants::size::kMaxNdetectors; ++iDet) {
+ aStream << std::setw(2) << std::setfill(' ') << this->GetNstationsGeometry(static_cast<L1DetectorID>(iDet)) << ' ';
}
- aStream << " | total = " << std::setw(2) << std::setfill(' ') << fNstationsGeometryTotal;
+ aStream << " | total = " << std::setw(2) << std::setfill(' ') << this->GetNstationsGeometry();
aStream << '\n';
aStream << indent << indentChar << indentChar << "Number of stations (Active): ";
- for (int idx = 0; idx < int(fNstationsActive.size()); ++idx) {
- aStream << std::setw(2) << std::setfill(' ') << fNstationsActive[idx] << ' ';
+ for (int iDet = 0; iDet < L1Constants::size::kMaxNdetectors; ++iDet) {
+ aStream << std::setw(2) << std::setfill(' ') << this->GetNstationsActive(static_cast<L1DetectorID>(iDet)) << ' ';
}
- aStream << " | total = " << std::setw(2) << std::setfill(' ') << fNstationsActiveTotal;
+ aStream << " | total = " << std::setw(2) << std::setfill(' ') << this->GetNstationsActive();
aStream << '\n';
aStream << indent << indentChar << indentChar << "Active station indexes: ";
for (int idx = 0; idx < fNstationsActiveTotal; ++idx) {
- aStream << std::setw(3) << std::setfill(' ') << fActiveStationGlobalIDs[idx] << ' ';
+ aStream << std::setw(3) << std::setfill(' ') << fvGeoToActiveMap[idx] << ' ';
}
aStream << '\n';
@@ -281,8 +272,7 @@ std::string L1Parameters::ToString(int verbosity, int indentLevel) const
aStream << indent << indentChar << "Triplets matching via MC: " << fDevIsMatchTripletsViaMc << '\n';
aStream << indent << indentChar << "Extend tracks with MC matching: " << fDevIsExtendTracksViaMc << '\n';
aStream << indent << indentChar << "Overlap hits matching via MC: " << fDevIsSuppressOverlapHitsViaMc << '\n';
-
- aStream << indent << indentChar << "Use hit search windows: " << fDevIsParSearchWUsed << '\n';
+ aStream << indent << indentChar << "Use hit search windows: " << fDevIsParSearchWUsed << '\n';
if (fDevIsParSearchWUsed) {
aStream << indent << "SEARCH WINDOWS:\n";
diff --git a/reco/L1/L1Algo/L1Parameters.h b/reco/L1/L1Algo/L1Parameters.h
index 03c235f250053cd887c93cc1800ba6f1932fb5ab..bc53ac71010d23a51002e341879bf1c318d26bd6 100644
--- a/reco/L1/L1Algo/L1Parameters.h
+++ b/reco/L1/L1Algo/L1Parameters.h
@@ -79,13 +79,6 @@ public:
/// 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; }
@@ -96,57 +89,52 @@ public:
int GetNstationsActive() const { return fNstationsActiveTotal; }
/// Gets number of active stations for given detector ID
- int GetNstationsActive(L1DetectorID detectorID) const
- {
- return fNstationsActive[static_cast<L1DetectorID_t>(detectorID)];
- }
+ int GetNstationsActive(L1DetectorID detectorID) const;
- /// Gets total number of stations, provided by setup geometry
- int GetNstationsGeometry() const { return fNstationsGeometryTotal; }
+ /// @brief Gets total number of stations, provided by setup geometry
+ int GetNstationsGeometry() const { return fvFirstGeoId.back(); }
- /// Gets number of stations, provided by setup geometry for given detector ID
+ /// @brief Gets number of stations, provided by setup geometry for given detector ID
int GetNstationsGeometry(L1DetectorID detectorID) const
{
- return fNstationsGeometry[static_cast<L1DetectorID_t>(detectorID)];
+ return fvFirstGeoId[static_cast<int>(detectorID) + 1] - fvFirstGeoId[static_cast<int>(detectorID)];
}
- /// Calculates global index of station among geometry (accounts for inactive stations)
- /// \param localIndex index of the detector subsystem module/station/layer provided by detector subsystem experts
- /// \param detectorID ID of the detector subsystem
+ /// @brief Calculates global index of station among geometry (accounts for inactive stations)
+ /// @param localIndex index of the detector subsystem module/station/layer provided by detector subsystem experts
+ /// @param detectorID ID of the detector subsystem
[[gnu::always_inline]] int GetStationIndexGeometry(int localIndex, L1DetectorID detectorID) const
{
- return localIndex
- + std::accumulate(fNstationsGeometry.cbegin(), fNstationsGeometry.cbegin() + static_cast<int>(detectorID),
- 0);
+ return fvLocalToGeoIdMap[fvFirstGeoId[static_cast<int>(detectorID)] + localIndex];
}
- /// Calculates global index of station used by track finder
- /// \param localIndex index of the detector subsystem module/station/layer provided by detector subsystem experts
- /// \param detectorID ID of the detector subsystem
+ /// @brief Calculates global index of station used by track finder
+ /// @param localIndex index of the detector subsystem module/station/layer provided by detector subsystem experts
+ /// @param detectorID ID of the detector subsystem
[[gnu::always_inline]] int GetStationIndexActive(int localIndex, L1DetectorID detectorID) const
{
- return fActiveStationGlobalIDs[GetStationIndexGeometry(localIndex, detectorID)];
+ return fvGeoToActiveMap[GetStationIndexGeometry(localIndex, detectorID)];
}
- /// Provides access to L1CAIteration vector (const)
+ /// @brief Provides access to L1CAIteration vector (const)
const L1IterationsContainer_t& GetCAIterations() const { return fCAIterations; }
- /// Provides number of iterations
+ /// @brief Provides number of iterations
int GetNcaIterations() const { return fCAIterations.size(); }
- /// Provides access to L1Stations container (const)
+ /// @brief Provides access to L1Stations container (const)
const L1StationsContainer_t& GetStations() const { return fStations; }
- /// Gets reference to the particular station
- /// \param iStation Index of station in the active stations container
+ /// @brief Gets reference to the particular station
+ /// @param iStation Index of station in the active stations container
const L1Station& GetStation(int iStation) const { return fStations[iStation]; }
- /// Gets a search window for a selected station and track group
- /// \note For a particular track finder iteration one can select a track group, which is defined by the minimal
- /// momentum of tracks (fast tracks, all tracks), their vertex (primary or secondary tracks), or by particle
- /// (electrons, muons, hadrons, etc.)
- /// \param iStation Global index of active station
- /// \param iTrackGr Index of a track group
+ /// @brief Gets a search window for a selected station and track group
+ /// @note For a particular track finder iteration one can select a track group, which is defined by the minimal
+ /// momentum of tracks (fast tracks, all tracks), their vertex (primary or secondary tracks), or by particle
+ /// (electrons, muons, hadrons, etc.)
+ /// @param iStation Global index of active station
+ /// @param iTrackGr Index of a track group
const L1SearchWindow& GetSearchWindow(int iStation, int iTrackGr) const
{
assert(iStation < fNstationsActiveTotal && iStation > 0);
@@ -154,13 +142,13 @@ public:
return fSearchWindows[iTrackGr * L1Constants::size::kMaxNstations + iStation];
}
- /// Gets reference to the array of station thickness map
+ /// @brief Gets reference to the array of station thickness map
const L1MaterialContainer_t& GetThicknessMaps() const { return fThickMap; }
- /// Gets material thickness in units of radiation length in a point on the XY plane for a selected station
- /// \param iStActive Global index of an active station
- /// \param xPos Position of the point in X dimension [cm]
- /// \param yPos Position of the point in Y dimension [cm]
+ /// @brief Gets material thickness in units of radiation length in a point on the XY plane for a selected station
+ /// @param iStActive Global index of an active station
+ /// @param xPos Position of the point in X dimension [cm]
+ /// @param yPos Position of the point in Y dimension [cm]
float GetMaterialThicknessScal(int iStActive, float xPos, float yPos) const
{
return fThickMap[iStActive].GetRadThickScal(xPos, yPos);
@@ -255,27 +243,37 @@ private:
/// Array of station thickness map
alignas(L1Constants::misc::kAlignment) L1MaterialContainer_t fThickMap {};
- int fNstationsGeometryTotal = -1; ///< total number of stations, provided by geometry
- int fNstationsActiveTotal = -1; ///< total number of active tracking stations
+ // ** Station layout arrays **
+ /// @brief First index of the station on the particular detector
+ /// The last element of the array corresponds to the total number of geometry stations
+ alignas(L1Constants::misc::kAlignment) std::array<int, L1Constants::size::kMaxNdetectors + 1> fvFirstGeoId {};
+
+ /// @brief Map of (local, det) to geo indices
+ ///
+ /// @note Usage:
+ /// iStGeo = fvGeoToLocalIdMap[fvFirstGeoId[iDet] + iStLocal];
+ /// geo index.
+ alignas(L1Constants::misc::kAlignment) std::array<int, L1Constants::size::kMaxNstations> fvLocalToGeoIdMap {};
+
+ /// @brief Map of geo to active indices
+ ///
+ /// The vector maps actual station index (which is defined by ) to the index of station in tracking. If the station
+ /// is inactive, its index is equal to -1.
+ /// @example Let stations 1 and 4 be inactive. Then:
+ /// geometry index: 0 1 2 3 4 5 6 7 8 9 0 0 0 0
+ /// active index: 0 -1 1 2 -1 3 4 5 6 7 0 0 0 0
+ alignas(L1Constants::misc::kAlignment) std::array<int, L1Constants::size::kMaxNstations> fvGeoToActiveMap {};
- /// Actual numbers of stations, provided by geometry. Index of an array element (except the last one) corresponds to a given
- /// L1DetectorID of the detector subsystem. The last array element corresponds to the total number of stations.
- alignas(L1Constants::misc::kAlignment) std::array<int, L1Constants::size::kMaxNdetectors> fNstationsGeometry = {};
+ /// @brief Map of active to geo indices
+ alignas(L1Constants::misc::kAlignment) std::array<int, L1Constants::size::kMaxNstations> fvActiveToGeoMap {};
- /// Numbers of stations, which are active in tracking. Index of an array element (except the last one) corresponds to a given
- /// L1DetectorID of the detector subsystem. The last array element corresponds to the total number of stations.
- alignas(L1Constants::misc::kAlignment) std::array<int, L1Constants::size::kMaxNdetectors> fNstationsActive = {};
+ alignas(L1Constants::misc::kAlignment) int fNstationsActiveTotal = -1; ///< total number of active tracking stations
- /// Map of the actual detector indexes to the active detector indexes
- /// The vector maps actual station index (which is defined by ) to the index of station in tracking. If the station is inactive,
- /// its index is equal to -1. Example: let stations 1 and 4 be inactive. Then:
- /// actual index: 0 1 2 3 4 5 6 7 8 9 0 0 0 0
- /// active index: 0 -1 1 2 -1 3 4 5 6 7 0 0 0 0
- alignas(L1Constants::misc::kAlignment) std::array<int, L1Constants::size::kMaxNstations> fActiveStationGlobalIDs {};
- /// Map of search windows vs. active station global index and tracks group
- /// The tracks group can be defined by minimum momentum (fast/all tracks), origin (primary/secondary) and particle type
- /// (electron, muon, all). Other options also can be added
+ /// @brief Map of search windows vs. active station global index and tracks group
+ ///
+ /// The tracks group can be defined by minimum momentum (fast/all tracks), origin (primary/secondary) and particle
+ /// type (electron, muon, all). Other options also can be added
alignas(L1Constants::misc::kAlignment)
std::array<L1SearchWindow, L1Constants::size::kMaxNstations* L1Constants::size::kMaxNtrackGroups> fSearchWindows =
{};
@@ -317,11 +315,10 @@ private:
ar& fStations;
ar& fThickMap;
- ar& fNstationsGeometryTotal;
+ ar& fvFirstGeoId;
+ ar& fvLocalToGeoIdMap;
+ ar& fvGeoToActiveMap;
ar& fNstationsActiveTotal;
- ar& fNstationsGeometry;
- ar& fNstationsActive;
- ar& fActiveStationGlobalIDs;
ar& fSearchWindows;
ar& fTrackingLevel;
diff --git a/reco/L1/qa/CbmCaInputQaMuch.cxx b/reco/L1/qa/CbmCaInputQaMuch.cxx
index bef4e3b95cf967fe64c392a08583d696a252f6c7..4cc613df15e3d47125746511c654a6ff1fb90983 100644
--- a/reco/L1/qa/CbmCaInputQaMuch.cxx
+++ b/reco/L1/qa/CbmCaInputQaMuch.cxx
@@ -366,7 +366,7 @@ void CbmCaInputQaMuch::FillHistograms()
fvph_n_points_per_hit[iSt]->Fill(nMCpoints);
fvph_n_points_per_hit[nSt]->Fill(nMCpoints);
- if (nMCpoints != 1) { continue; } // ??
+ if (nMCpoints != 1) { continue; }
// The best link to in the match (probably, the cut on nMCpoints is meaningless)
const auto& bestPointLink = pHitMatch->GetMatchedLink();
@@ -481,7 +481,7 @@ void CbmCaInputQaMuch::FillHistograms()
int iFile = fpMCEventList->GetFileIdByIndex(iE);
int iEvent = fpMCEventList->GetEventIdByIndex(iE);
int nPoints = fpMCPoints->Size(iFile, iEvent);
-
+ std::cout << "points:\n";
for (int iP = 0; iP < nPoints; ++iP) {
const auto* pMCPoint = dynamic_cast<const CbmMuchPoint*>(fpMCPoints->Get(iFile, iEvent, iP));
LOG_IF(fatal, !pMCPoint) << fName << ": MC point does not exist for iFile = " << iFile
@@ -496,9 +496,7 @@ void CbmCaInputQaMuch::FillHistograms()
double yMC = pMCPoint->GetYIn();
double rMC = sqrt(xMC * xMC + yMC * yMC);
- // Conditions under which point is accounted as reconstructed: point
bool ifPointHasHits = (vNofHitsPerPoint[iE][iP] > 0);
-
fvpe_reco_eff_vs_xy[iSt]->Fill(ifPointHasHits, xMC, yMC);
fvpe_reco_eff_vs_xy[nSt]->Fill(ifPointHasHits, xMC, yMC);