-
Sergei Zharko authoredSergei Zharko authored
CaParameters.cxx 14.66 KiB
/* Copyright (C) 2022-2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
SPDX-License-Identifier: GPL-3.0-only
Authors: Sergey Gorbunov, Sergei Zharko [committer] */
/// \file L1Parameters.cxx
/// \brief Parameter container for the L1Algo library
/// \since 10.02.2022
/// \author S.Zharko <s.zharko@gsi.de>
#include "CaParameters.h"
#include <iomanip>
#include "AlgoFairloggerCompat.h"
using cbm::algo::ca::EDetectorID;
using cbm::algo::ca::Iteration;
using cbm::algo::ca::MaterialMap;
using cbm::algo::ca::Parameters;
using cbm::algo::ca::utils::CheckSimdVectorEquality;
// ---------------------------------------------------------------------------------------------------------------------
//
Parameters::Parameters()
{
fvGeoToActiveMap.fill(-1); // by default, all stations are inactive, thus all the IDs must be -1
}
// ---------------------------------------------------------------------------------------------------------------------
//
void Parameters::CheckConsistency() const
{
LOG(info) << "Consistency test for L1 parameters object... ";
/*
* Array of active station IDs
*
* In the array of active station IDs,
* (i) sum of elements, which are not equal -1, must be equal the number of stations,
* (ii) subsequence of all the elements, which are not equal -1, must be a gapless subset of integer numbers starting with 0
*/
auto filterInactiveStationIDs = [](int x) { return x != -1; };
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 "
<< "with the total number of stations (" << nStationsCheck << " vs. " << fNstationsActiveTotal << ' '
<< "expected)";
throw std::logic_error(msg.str());
}
// Check, if the subsequence of all the elements, which are not equal -1, must be a gapless subset of integer numbers
// starting from 0. If it is, the testValue in the end must be equal the number of non -1 elements
std::vector<int> idsCheck(
nStationsCheck); // temporary vector containing a sequence of active station IDs without -1 elements
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;
}
if (!isStationIDsOk) {
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 : fvGeoToActiveMap) {
msg << std::setw(3) << std::setfill(' ') << id << ' ';
}
throw std::logic_error(msg.str());
}
/*
* Check magnetic field flags of the stations
*
* In a current version of tracking there are three configurations possible to be proceeded:
* A. All the stations are inside magnetic field
* B. There is no magnetic field in a setup
* C. All the first stations are inside magnetic field, all the last stations are outside the field
* In all the cases the fieldStatus flags should be sorted containing all non-zero elements in the beginning
* (representing stations placed into magnetic field) and all zero elements in the end of z-axis.
*/
bool ifFieldStatusFlagsOk = std::is_sorted(
fStations.cbegin(), fStations.cbegin() + fNstationsActiveTotal,
[&](const ca::Station& lhs, const ca::Station& rhs) { return bool(lhs.fieldStatus) > bool(rhs.fieldStatus); });
if (!ifFieldStatusFlagsOk) {
std::stringstream msg;
msg << "L1Parameters: invalid object condition: L1 tracking is impossible for a given field configuration:\n";
for (int iSt = 0; iSt < fNstationsActiveTotal; ++iSt) {
msg << "- station ID: " << iSt << ", field status: " << fStations[iSt].fieldStatus << '\n';
}
throw std::logic_error(msg.str());
}
/*
* Check target position SIMD vector
*/
utils::CheckSimdVectorEquality(fTargetPos[0], "L1Parameters: target position x");
utils::CheckSimdVectorEquality(fTargetPos[1], "L1Parameters: target position y");
utils::CheckSimdVectorEquality(fTargetPos[2], "L1Parameters: target position z");
/*
* Check vertex field region and value objects at primary vertex
*/
fVertexFieldValue.CheckConsistency();
fVertexFieldRegion.CheckConsistency();
/*
* Check if each station object is consistent itself, and if all of them are placed after the target
* NOTE: If a station was not set up, it is accounted inconsistent (uninitialized). In the array of stations there are uninitialized
* stations possible (with id > NstationsActiveTotal), thus one should NOT run the loop above over all the stations in array
* but only until *(fNstationsActive.cend() - 1) (== NstationsActiveTotal).
* TODO: Probably, we should introduce methods, which check the consistency of fully initialized objects such as ca::Station,
* L1MaterialInfo, etc. (S.Zharko)
*/
for (int iSt = 0; iSt < fNstationsActiveTotal; ++iSt) {
fStations[iSt].CheckConsistency();
if (fStations[iSt].fZ[0] < fTargetPos[2][0]) {
std::stringstream msg;
msg << "L1Parameters: station with global ID = " << iSt << " is placed before target "
<< "(z_st = " << fStations[iSt].fZ[0] << " [cm] < z_targ = " << fTargetPos[2][0] << " [cm])";
throw std::logic_error(msg.str());
}
}
/*
* Check thickness maps
* NOTE: If a ca::MaterialMap map was not set up, it is accounted inconsistent (uninitialized). In the array of thickness maps for each
* there are uninitialized elements possible (with id > NstationsActiveTotal), thus one should NOT run the loop above over
* all the stations in array but only until *(fNstationsActive.cend() - 1) (== NstationsActiveTotal).
*/
for (int iSt = 0; iSt < fNstationsActiveTotal; ++iSt) {
fThickMap[iSt].CheckConsistency();
}
/*
* Check iterations sequence
* 1. Number of iterations should be larger then zero
* 2. Each iteration should contain values within predefined limits
* 3. Number of iterations with TrackFromTriplets flag turned on no more then 1
* 4. If the TrackFromTriplets iteration exists, it should be the last one in the sequence
*/
{
int nIterations = fCAIterations.size();
if (!nIterations) {
std::stringstream msg;
msg << "L1Parameters: 0 track finder iterations were found. Please, define at least one iteration";
throw std::logic_error(msg.str());
}
std::string names = "";
for (const auto& iter : fCAIterations) {
if (!iter.Check()) { names += iter.GetName() + " "; }
}
if (names.size()) {
std::stringstream msg;
msg << "L1Parameters: some parameters are out of range for the following iterations: " << names;
throw std::logic_error(msg.str());
}
nIterations = std::count_if(fCAIterations.begin(), fCAIterations.end(),
[=](const Iteration& it) { return it.GetTrackFromTripletsFlag(); });
if (nIterations > 1) {
std::stringstream msg;
msg << "L1Parameters: found " << nIterations << " iterations with GetTrackFromTripletsFlag() == true:\n";
for (const auto& iter : fCAIterations) {
if (iter.GetTrackFromTripletsFlag()) { msg << '\t' << "- " << iter.GetName() << '\n'; }
}
msg << "Only the one iteration can have GetTrackFromTripletsFlag() == true, and this iteration should be ";
msg << "the last one";
throw std::logic_error(msg.str());
}
if (nIterations == 1 && !(fCAIterations.end() - 1)->GetTrackFromTripletsFlag()) {
std::stringstream msg;
msg << "L1Parameters: iteration with GetTrackFromTripletsFlag() == true is not the last in a sequence. ";
msg << "The GetTrackFromTripletsFlag() value in the iterations sequence: \n";
for (const auto& iter : fCAIterations) {
msg << '\t' << "- " << std::setw(15) << std::setfill(' ') << iter.GetName() << ' ';
msg << std::setw(6) << std::setfill(' ') << iter.GetTrackFromTripletsFlag() << '\n';
}
throw std::logic_error(msg.str());
}
}
LOG(info) << "Consistency test for L1 parameters object... \033[1;32mpassed\033[0m";
}
// ---------------------------------------------------------------------------------------------------------------------
//
int Parameters::GetNstationsActive(EDetectorID 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 Parameters::Print(int /*verbosityLevel*/) const { LOG(info) << ToString(); }
// ---------------------------------------------------------------------------------------------------------------------
//
std::string Parameters::ToString(int verbosity, int indentLevel) const
{
using namespace constants;
using std::setfill;
using std::setw;
std::stringstream msg {};
if (verbosity < 1) { return msg.str(); }
constexpr char indentCh = '\t';
std::string indent(indentLevel, indentCh);
/// Sets red/green color dependently from flag
auto ClrFlg = [&](bool flag) { return (flag ? clrs::GN : clrs::RD); };
msg << " ----- CA parameters list -----\n";
msg << indent << clrs::CLb << "COMPILE TIME CONSTANTS:\n" << clrs::CL;
msg << indent << indentCh << "Bits to code one station: " << constants::size::StationBits << '\n';
msg << indent << indentCh << "Bits to code one triplet: " << constants::size::TripletBits << '\n';
msg << indent << indentCh << "Max number of detectors: " << constants::size::MaxNdetectors << '\n';
msg << indent << indentCh << "Max number of stations: " << constants::size::MaxNstations << '\n';
msg << indent << indentCh << "Max number of triplets: " << constants::size::MaxNtriplets << '\n';
msg << indent << clrs::CLb << "RUNTIME CONSTANTS:\n" << clrs::CL;
msg << indent << indentCh << "Random seed: " << fRandomSeed << '\n';
msg << indent << indentCh << "Max number of doublets per singlet: " << fMaxDoubletsPerSinglet << '\n';
msg << indent << indentCh << "Max number of triplets per doublet: " << fMaxTripletPerDoublets << '\n';
msg << indent << indentCh << "Ghost suppresion: " << fGhostSuppression << '\n';
msg << indent << clrs::CLb << "CA TRACK FINDER ITERATIONS:\n" << clrs::CL;
msg << Iteration::ToTableFromVector(fCAIterations);
msg << indent << clrs::CLb << "GEOMETRY:\n" << clrs::CL;
msg << indent << indentCh << clrs::CLb << "TARGET:\n" << clrs::CL;
msg << indent << indentCh << indentCh << "Position:\n";
for (int dim = 0; dim < 3 /*nDimensions*/; ++dim) {
msg << indent << indentCh << indentCh << indentCh << char(120 + dim) << " = " << fTargetPos[dim][0] << " cm\n";
}
msg << indent << indentCh << clrs::CLb << "NUMBER OF STATIONS:\n" << clrs::CL;
msg << indent << indentCh << indentCh << "Number of stations (Geometry): ";
for (int iDet = 0; iDet < constants::size::MaxNdetectors; ++iDet) {
msg << setw(2) << setfill(' ') << this->GetNstationsGeometry(static_cast<EDetectorID>(iDet)) << ' ';
}
msg << " | total = " << setw(2) << setfill(' ') << this->GetNstationsGeometry();
msg << '\n';
msg << indent << indentCh << indentCh << "Number of stations (Active): ";
for (int iDet = 0; iDet < constants::size::MaxNdetectors; ++iDet) {
msg << setw(2) << setfill(' ') << this->GetNstationsActive(static_cast<EDetectorID>(iDet)) << ' ';
}
msg << " | total = " << setw(2) << setfill(' ') << this->GetNstationsActive();
msg << '\n' << indent << indentCh << indentCh << clrs::CL << "Geometry station indices: ";
for (int iStGeo = 0; iStGeo < this->GetNstationsGeometry(); ++iStGeo) {
bool isActive = fvGeoToActiveMap[iStGeo] != -1;
msg << ClrFlg(isActive) << setw(3) << setfill(' ') << iStGeo << ' ';
}
msg << '\n' << indent << indentCh << indentCh << clrs::CL << "Local station indices: ";
for (int iStGeo = 0; iStGeo < this->GetNstationsGeometry(); ++iStGeo) {
bool isActive = fvGeoToActiveMap[iStGeo] != -1;
msg << ClrFlg(isActive) << setw(3) << setfill(' ') << static_cast<int>(fvGeoToLocalIdMap[iStGeo].second) << ' ';
}
msg << '\n' << indent << indentCh << indentCh << clrs::CL << "Detector indices: ";
for (int iStGeo = 0; iStGeo < this->GetNstationsGeometry(); ++iStGeo) {
bool isActive = fvGeoToActiveMap[iStGeo] != -1;
msg << ClrFlg(isActive) << setw(3) << setfill(' ') << static_cast<int>(fvGeoToLocalIdMap[iStGeo].first) << ' ';
}
msg << '\n' << indent << indentCh << indentCh << clrs::CL << "Active station indices: ";
for (int iStGeo = 0; iStGeo < this->GetNstationsGeometry(); ++iStGeo) {
bool isActive = fvGeoToActiveMap[iStGeo] != -1;
msg << ClrFlg(isActive) << setw(3) << setfill(' ') << fvGeoToActiveMap[iStGeo] << ' ';
}
msg << clrs::CL << '\n';
msg << indent << indentCh << clrs::CLb << "STATIONS:\n" << clrs::CL;
msg << indent << indentCh << setw(9) << setfill(' ') << "Active ID" << ' ';
msg << fStations[0].ToString(1, 0, true) << '\n';
for (int iStAct = 0; iStAct < this->GetNstationsActive(); ++iStAct) {
msg << indent << indentCh << setw(9) << setfill(' ') << iStAct << ' ';
msg << fStations[iStAct].ToString(verbosity, 0) << '\n';
}
msg << indent << clrs::CLb << "DEV FLAGS:" << clrs::CL << " (for debug only)\n";
msg << indent << indentCh << "Hits search area is ignored: " << fDevIsIgnoreHitSearchAreas << '\n';
msg << indent << indentCh << "Non-approx. field is used: " << fDevIsMatchDoubletsViaMc << '\n';
msg << indent << indentCh << "Doublets matching via MC: " << fDevIsMatchDoubletsViaMc << '\n';
msg << indent << indentCh << "Triplets matching via MC: " << fDevIsMatchTripletsViaMc << '\n';
msg << indent << indentCh << "Extend tracks with MC matching: " << fDevIsExtendTracksViaMc << '\n';
msg << indent << indentCh << "Overlap hits matching via MC: " << fDevIsSuppressOverlapHitsViaMc << '\n';
msg << indent << indentCh << "Use hit search windows: " << fDevIsParSearchWUsed << '\n';
if (fDevIsParSearchWUsed) {
msg << indent << "SEARCH WINDOWS:\n";
for (int iSt = 1; iSt < fNstationsActiveTotal; ++iSt) {
for (int iIter = 0; iIter < (int) fCAIterations.size(); ++iIter) {
msg << indent << "- station: " << iSt << ", iteration: " << iIter << '\n';
msg << GetSearchWindow(iSt, iIter).ToString() << '\n';
}
}
}
msg << '\n';
return msg.str();
}