diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx
index a7cbab2f40fa4c189c1dc2302bc5d23aaa689da2..777a0c63e2417c2de37813f7c40729c80b4771c7 100644
--- a/reco/L1/CbmL1.cxx
+++ b/reco/L1/CbmL1.cxx
@@ -547,11 +547,40 @@ InitStatus CbmL1::Init()
// LOG(error) << "-E- CbmL1: Read geometry from file " << fSTAPDataDir + "geo_algo.txt was NOT successful.";
}
+
+ /****************************
+ ** Material budget input **
+ ****************************/
+
+ // NOTE: std::vector of material tables. Vector sizes correspond to number of stations provided by geometry, i.e. both
+ // active and inactive station are represented in the folloving vectors
+ auto materialTableMvd = ReadMaterialBudget(L1DetectorID::kMvd);
+ auto materialTableSts = ReadMaterialBudget(L1DetectorID::kSts);
+ auto materialTableMuch = ReadMaterialBudget(L1DetectorID::kMuch);
+ auto materialTableTrd = ReadMaterialBudget(L1DetectorID::kTrd);
+ auto materialTableTof = ReadMaterialBudget(L1DetectorID::kTof);
+
+ /* User corrections (optional) */
+ auto correctionMvd = [this](L1Material& material, const L1MaterialInfo& homogenious) {
+ this->ApplyCorrectionToMaterialMap<L1DetectorID::kMvd>(material, homogenious);
+ };
+ auto correctionSts = [this](L1Material& material, const L1MaterialInfo& homogenious) {
+ this->ApplyCorrectionToMaterialMap<L1DetectorID::kSts>(material, homogenious);
+ };
+ auto correctionMuch = [this](L1Material& material, const L1MaterialInfo& homogenious) {
+ this->ApplyCorrectionToMaterialMap<L1DetectorID::kMuch>(material, homogenious);
+ };
+ auto correctionTrd = [this](L1Material& material, const L1MaterialInfo& homogenious) {
+ this->ApplyCorrectionToMaterialMap<L1DetectorID::kTrd>(material, homogenious);
+ };
+ auto correctionTof = [this](L1Material& material, const L1MaterialInfo& homogenious) {
+ this->ApplyCorrectionToMaterialMap<L1DetectorID::kTof>(material, homogenious);
+ };
+
/***************************************
** Stations geometry initialization **
***************************************/
-
/*** MVD stations info ***/
for (int iSt = 0; iSt < NMvdStationsGeom; ++iSt) { // NOTE: example using in-stack defined objects
CbmMvdDetector* mvdDetector = CbmMvdDetector::Instance();
@@ -567,7 +596,8 @@ InitStatus CbmL1::Init()
stationInfo.SetZ(t.z);
auto thickness = t.dz;
auto radLength = t.RadLength;
- stationInfo.SetMaterial(thickness, radLength);
+ stationInfo.SetMaterialSimple(thickness, radLength);
+ stationInfo.SetMaterialMap(std::move(materialTableMvd[iSt]), correctionMvd);
stationInfo.SetXmax(t.R);
stationInfo.SetYmax(t.R);
stationInfo.SetRmin(t.r);
@@ -600,7 +630,8 @@ InitStatus CbmL1::Init()
stationInfo.SetRmax(stsXmax > stsYmax ? stsXmax : stsYmax);
auto thickness = cbmSts->GetSensorD();
auto radLength = cbmSts->GetRadLength();
- stationInfo.SetMaterial(thickness, radLength);
+ stationInfo.SetMaterialSimple(thickness, radLength);
+ stationInfo.SetMaterialMap(std::move(materialTableSts[iSt]), correctionSts);
// Setup strips geometry
fscal stsFrontPhi = cbmSts->GetSensorRotation() + cbmSts->GetSensorStereoAngle(0) * TMath::Pi() / 180.;
fscal stsBackPhi = cbmSts->GetSensorRotation() + cbmSts->GetSensorStereoAngle(1) * TMath::Pi() / 180.;
@@ -627,7 +658,8 @@ InitStatus CbmL1::Init()
stationInfo.SetZ(layer->GetZ());
auto thickness = layer->GetDz();
auto radLength = 0.; // Why 0??? (S.Zharko)
- stationInfo.SetMaterial(thickness, radLength);
+ stationInfo.SetMaterialSimple(thickness, radLength);
+ stationInfo.SetMaterialMap(std::move(materialTableMuch[iSt]), correctionMuch);
stationInfo.SetXmax(100.);
stationInfo.SetYmax(100.);
stationInfo.SetRmin(10.);
@@ -647,7 +679,7 @@ InitStatus CbmL1::Init()
int trdModuleID = fTrdDigiPar->GetModuleId(iSt);
CbmTrdParModDigi* module = (CbmTrdParModDigi*) fTrdDigiPar->GetModulePar(trdModuleID);
auto stationInfo = L1BaseStationInfo(L1DetectorID::kTrd, iSt);
- int stationType = (iSt == 2 || iSt == 4) ? 6 : 3; // Is used somewhere??
+ int stationType = (iSt == 1 || iSt == 3) ? 6 : 3; // Is used somewhere??
stationInfo.SetStationType(stationType);
stationInfo.SetTimeInfo(1);
stationInfo.SetTimeResolution(10.);
@@ -655,7 +687,8 @@ InitStatus CbmL1::Init()
stationInfo.SetZ(module->GetZ());
auto thickness = 2. * module->GetSizeZ();
auto radLength = 1.6;
- stationInfo.SetMaterial(thickness, radLength);
+ stationInfo.SetMaterialSimple(thickness, radLength);
+ stationInfo.SetMaterialMap(std::move(materialTableTrd[iSt]), correctionTrd);
stationInfo.SetXmax(module->GetSizeX());
stationInfo.SetYmax(module->GetSizeY());
stationInfo.SetRmin(0.);
@@ -683,7 +716,8 @@ InitStatus CbmL1::Init()
stationInfo.SetZ(TofStationZ[iSt]);
auto thickness = 10.;
auto radLength = 2.;
- stationInfo.SetMaterial(thickness, radLength);
+ stationInfo.SetMaterialSimple(thickness, radLength);
+ stationInfo.SetMaterialMap(std::move(materialTableTof[iSt]), correctionTof);
stationInfo.SetXmax(20.);
stationInfo.SetYmax(20.);
stationInfo.SetRmin(0.);
@@ -892,19 +926,6 @@ InitStatus CbmL1::Init()
algo->Init(fUseHitErrors, fTrackingMode, fMissingHits);
- /*********************************
- ** Material map initialization **
- *********************************/
-
- // TODO: Move it to the L1InitManager (S.Zharko, 15.05.2022)
- algo->fRadThick.reset(algo->GetNstations());
-
- if (fUseMVD) { ReadMaterialTables(L1DetectorID::kMvd); }
- ReadMaterialTables(L1DetectorID::kSts);
- if (fUseMUCH) { ReadMaterialTables(L1DetectorID::kMuch); }
- if (fUseTRD) { ReadMaterialTables(L1DetectorID::kTrd); }
- if (fUseTOF) { ReadMaterialTables(L1DetectorID::kTof); }
-
return kSUCCESS;
}
@@ -2218,7 +2239,6 @@ void CbmL1::ReadMaterialTables(L1DetectorID detectorID)
int nBins = hStaRadLen->GetNbinsX();
float rMax = hStaRadLen->GetXaxis()->GetXmax();
algo->fRadThick[iStActive].SetBins(nBins, rMax);
- algo->fRadThick[iStActive].table.resize(nBins);
float hole = 0.f;
switch (detectorID) {
@@ -2231,16 +2251,16 @@ void CbmL1::ReadMaterialTables(L1DetectorID detectorID)
double averageRadThick = 0.;
double counter = 0.;
for (int iBinX = 0; iBinX < nBins; ++iBinX) {
- algo->fRadThick[iStActive].table[iBinX].resize(nBins);
if (detectorID == L1DetectorID::kTof) { hole = 0.0015f; } // TODO: Why? (S.Zharko)
for (int iBinY = 0; iBinY < nBins; ++iBinY) {
- algo->fRadThick[iStActive].table[iBinX][iBinY] = 0.01 * hStaRadLen->GetBinContent(iBinX, iBinY);
+ algo->fRadThick[iStActive].SetRadThick(iBinX, iBinY, 0.01 * hStaRadLen->GetBinContent(iBinX, iBinY));
/* Specific corrections for each detector type */
+ // TODO: Is it possible to account for these corrections in the moment of the maps creation? (S.Zharko)
switch (detectorID) {
case L1DetectorID::kMvd:
// Correction for holes in the material map
- if (algo->fRadThick[iStActive].table[iBinX][iBinY]
+ if (algo->fRadThick[iStActive].GetRadThick(iBinX, iBinY)
< algo->GetStations()[iStActive].materialInfo.RadThick[0]) {
if (iBinY > 0 && iBinY < nBins - 1) {
algo->fRadThick[iStActive].table[iBinX][iBinY] =
@@ -2255,10 +2275,10 @@ void CbmL1::ReadMaterialTables(L1DetectorID detectorID)
}
break;
case L1DetectorID::kSts:
- if (algo->fRadThick[iStActive].table[iBinX][iBinY]
+ if (algo->fRadThick[iStActive].GetRadThick(iBinX, iBinY)
< algo->GetStations()[iStActive].materialInfo.RadThick[0]) {
- algo->fRadThick[iStActive].table[iBinX][iBinY] =
- algo->GetStations()[iStActive].materialInfo.RadThick[0];
+ algo->fRadThick[iStActive].SetRadThick(iBinX, iBinY,
+ algo->GetStations()[iStActive].materialInfo.RadThick[0]);
}
break;
case L1DetectorID::kMuch:
@@ -2266,6 +2286,7 @@ void CbmL1::ReadMaterialTables(L1DetectorID detectorID)
case L1DetectorID::kTof:
// Correction for holes in the material map
if (iBinY > 0 && iBinY < nBins - 1) {
+<<<<<<< HEAD
algo->fRadThick[iStActive].table[iBinX][iBinY] =
0.01
* TMath::Min(hStaRadLen->GetBinContent(iBinX, iBinY - 1),
@@ -2279,8 +2300,8 @@ void CbmL1::ReadMaterialTables(L1DetectorID detectorID)
algo->fRadThick[iStActive].table[iBinX][iBinY] = hole;
}
break;
- } // switch (detectorID)
- averageRadThick += algo->fRadThick[iStActive].table[iBinX][iBinY];
+ } // switoch (detectorID)
+ averageRadThick += algo->fRadThick[iStActive].GetRadThick(iBinX, iBinY);
++counter;
} // iBinY
} // iBinX
@@ -2300,11 +2321,54 @@ void CbmL1::ReadMaterialTables(L1DetectorID detectorID)
int iStActive = fpInitManager->GetStationIndexActive(iSt, detectorID);
if (iStActive == -1) { continue; }
algo->fRadThick[iStActive].SetBins(1, 100);
- algo->fRadThick[iStActive].table.resize(1);
- algo->fRadThick[iStActive].table[0].resize(1);
- algo->fRadThick[iStActive].table[0][0] = algo->GetStations()[iStActive].materialInfo.RadThick[0];
+ algo->fRadThick[iStActive].SetRadThick(0, 0, algo->GetStations()[iStActive].materialInfo.RadThick[0]);
LOG(info) << " Material for " << GetDetectorName(detectorID) << ": "
<< algo->GetStations()[iStActive].materialInfo.RadThick[0];
} // iSt
}
}
+
+std::vector<L1Material> CbmL1::ReadMaterialBudget(L1DetectorID detectorID)
+{
+ std::vector<L1Material> result {};
+ if (fMatBudgetFileName.find(detectorID) != fMatBudgetFileName.end()) {
+ auto* oldFile = gFile;
+ auto* oldDir = gDirectory;
+
+ auto rlFile = TFile(fMatBudgetFileName.at(detectorID));
+ if (rlFile.IsZombie()) { LOG(fatal) << "File " << fMatBudgetFileName.at(detectorID) << " is zombie!"; }
+ else {
+ LOG(info) << "Reading material budget for " << GetDetectorName(detectorID) << " from file "
+ << fMatBudgetFileName.at(detectorID);
+ }
+
+ result.resize(fpInitManager->GetNstationsGeometry(detectorID));
+ TString stationNamePrefix = "Radiation Thickness [%], Station";
+
+ // NOTE: Loop over geometry stations. We probably do not know which stations are active/inactive (S.Zharko)
+ for (int iSt = 0; iSt < fpInitManager->GetNstationsGeometry(detectorID); ++iSt) {
+ // TODO: Unify material table names (S.Zharko)
+ TString stationName = stationNamePrefix + (detectorID == L1DetectorID::kMvd ? iSt : iSt + 1);
+ auto* hStaRadLen = dynamic_cast<TProfile2D*>(rlFile.Get(stationName));
+ if (!hStaRadLen) {
+ LOG(fatal) << "CbmL1: material budget profile " << stationName << " does not exist in file "
+ << fMatBudgetFileName.at(detectorID);
+ }
+ int nBins = hStaRadLen->GetNbinsX();
+ float rMax = hStaRadLen->GetXaxis()->GetXmax();
+ result[iSt].SetBins(nBins, rMax);
+ for (int iBinX = 0; iBinX < nBins; ++iBinX) {
+ for (int iBinY = 0; iBinY < nBins; ++iBinY) {
+ result[iSt].SetRadThick(iBinX, iBinY, 0.01 * hStaRadLen->GetBinContent(iBinX, iBinY));
+ } // iBinX
+ } // iBinY
+ LOG(info) << "- station " << iSt;
+ } // iSt
+ gFile = oldFile;
+ gDirectory = oldDir;
+ } // if mat budget file found
+ else {
+ LOG(warn) << "No material budget file is found for " << GetDetectorName(detectorID);
+ }
+ return result;
+}
diff --git a/reco/L1/CbmL1.h b/reco/L1/CbmL1.h
index dbd39db930f74f6e2c991f9fae60f33a9c55b247..917d7882dab58a0bc3a583db5f78025f48d2b507 100644
--- a/reco/L1/CbmL1.h
+++ b/reco/L1/CbmL1.h
@@ -200,6 +200,8 @@ public:
/// \param detectorID ID of a detector subsystem
void ReadMaterialTables(L1DetectorID detectorID);
+ std::vector<L1Material> ReadMaterialBudget(L1DetectorID detectorID);
+
void SetExtrapolateToTheEndOfSTS(bool b) { fExtrapolateToTheEndOfSTS = b; }
void SetLegacyEventMode(bool b) { fLegacyEventMode = b; }
diff --git a/reco/L1/L1Algo/L1Algo.cxx b/reco/L1/L1Algo/L1Algo.cxx
index 1d1d4f40d37517f65d94d47f85f1395592fd282d..70574c349dd1b5ce0719b4d1664ee66bdaa9492f 100644
--- a/reco/L1/L1Algo/L1Algo.cxx
+++ b/reco/L1/L1Algo/L1Algo.cxx
@@ -108,6 +108,8 @@ void L1Algo::Init(const bool UseHitErrors, const TrackingMode mode, const bool M
// Fill L1Station array
fInitManager.TransferL1StationArray(fStations);
+ fRadThick.reset(fNstations);
+ fInitManager.TransferL1MaterialArray(fRadThick);
fTrackingLevel = fInitManager.GetTrackingLevel();
fGhostSuppression = fInitManager.GetGhostSuppression();
diff --git a/reco/L1/L1Algo/L1Algo.h b/reco/L1/L1Algo/L1Algo.h
index ef4e05ebf2bb35177da315f18fdf8932143cfb8d..54e65b5db0b04994514f4570f8f4005a9bffc8a2 100644
--- a/reco/L1/L1Algo/L1Algo.h
+++ b/reco/L1/L1Algo/L1Algo.h
@@ -227,18 +227,40 @@ private:
int fNfieldStations {0}; ///< number of stations in the field region
alignas(16) L1StationsArray_t fStations; ///< array of L1Station objects
+ L1Vector<L1Material> fRadThick {"fRadThick"}; // material for each station
+
public:
/// Gets total number of stations used in tracking
int GetNstations() const { return fNstations; }
+
/// Gets number of stations before the pipe (MVD stations in CBM)
int GetNstationsBeforePipe() const { return fNstationsBeforePipe; }
+
/// Gets number of stations situated in field region (MVD + STS in CBM)
int GetNfieldStations() const { return fNfieldStations; }
+
/// Gets reference to the stations array
const L1StationsArray_t& GetStations() const { return fStations; }
+ /// 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 GetMaterialThickness(int iStActive, float xPos, float yPos) const
+ {
+ return fRadThick[iStActive].GetRadThick(xPos, yPos);
+ }
+
+ /// 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] (SIMDized vector)
+ /// \param yPos Position of the point in Y dimension [cm] (SIMDized vector)
+ fvec GetMaterialThickness(int iStActive, fvec xPos, fvec yPos) const
+ {
+ return fRadThick[iStActive].GetRadThick(xPos, yPos);
+ }
+
public:
- L1Vector<L1Material> fRadThick {"fRadThick"}; // material for each station
int NStsStrips {0}; ///> number of strips
L1Vector<L1Hit>* vStsHits {nullptr}; ///> hits as a combination of front-, backstrips and z-position
diff --git a/reco/L1/L1Algo/L1BaseStationInfo.cxx b/reco/L1/L1Algo/L1BaseStationInfo.cxx
index 0d12f7d31e45037a01699790aad803b3ab0971f9..e52b621c17cadc31a4cbed1f73fad6186265438f 100644
--- a/reco/L1/L1Algo/L1BaseStationInfo.cxx
+++ b/reco/L1/L1Algo/L1BaseStationInfo.cxx
@@ -30,14 +30,14 @@
// CONSTRUCTORS AND DESTRUCTOR
//
-//----------------------------------------------------------------------------------------------------------------------//
+//------------------------------------------------------------------------------------------------------------------------
//
L1BaseStationInfo::L1BaseStationInfo() noexcept
{
LOG(debug) << "L1BaseStationInfo: Default constructor called for " << this << '\n'; // Temporary
}
-//----------------------------------------------------------------------------------------------------------------------//
+//------------------------------------------------------------------------------------------------------------------------
//
L1BaseStationInfo::L1BaseStationInfo(L1DetectorID detectorID, int stationID) noexcept
: fDetectorID(detectorID)
@@ -48,14 +48,14 @@ L1BaseStationInfo::L1BaseStationInfo(L1DetectorID detectorID, int stationID) noe
fInitController.SetFlag(EInitKey::kStationID);
}
-//----------------------------------------------------------------------------------------------------------------------//
+//------------------------------------------------------------------------------------------------------------------------
//
L1BaseStationInfo::~L1BaseStationInfo() noexcept
{
LOG(debug) << "L1BaseStationInfo: Destructor called for " << this << '\n'; // Temporary
}
-//----------------------------------------------------------------------------------------------------------------------//
+//------------------------------------------------------------------------------------------------------------------------
//
L1BaseStationInfo::L1BaseStationInfo(const L1BaseStationInfo& other) noexcept
: fDetectorID(other.fDetectorID)
@@ -65,12 +65,13 @@ L1BaseStationInfo::L1BaseStationInfo(const L1BaseStationInfo& other) noexcept
, fYmax(other.fYmax)
, fZPos(other.fZPos)
, fL1Station(other.fL1Station)
+ , fThicknessMap(other.fThicknessMap)
, fInitController(other.fInitController)
{
LOG(debug) << "L1BaseStationInfo: Copy constructor called: " << &other << " was copied into " << this;
}
-//----------------------------------------------------------------------------------------------------------------------//
+//------------------------------------------------------------------------------------------------------------------------
//
L1BaseStationInfo::L1BaseStationInfo(L1BaseStationInfo&& other) noexcept
{
@@ -78,7 +79,7 @@ L1BaseStationInfo::L1BaseStationInfo(L1BaseStationInfo&& other) noexcept
this->Swap(other);
}
-//----------------------------------------------------------------------------------------------------------------------//
+//------------------------------------------------------------------------------------------------------------------------
//
L1BaseStationInfo& L1BaseStationInfo::operator=(const L1BaseStationInfo& other) noexcept
{
@@ -87,7 +88,7 @@ L1BaseStationInfo& L1BaseStationInfo::operator=(const L1BaseStationInfo& other)
return *this;
}
-//----------------------------------------------------------------------------------------------------------------------//
+//------------------------------------------------------------------------------------------------------------------------
//
L1BaseStationInfo& L1BaseStationInfo::operator=(L1BaseStationInfo&& other) noexcept
{
@@ -104,7 +105,7 @@ L1BaseStationInfo& L1BaseStationInfo::operator=(L1BaseStationInfo&& other) noexc
// BASIC METHODS
//
-//----------------------------------------------------------------------------------------------------------------------//
+//------------------------------------------------------------------------------------------------------------------------
//
void L1BaseStationInfo::Print(int verbosity) const
{
@@ -125,7 +126,7 @@ void L1BaseStationInfo::Print(int verbosity) const
}
}
-//----------------------------------------------------------------------------------------------------------------------//
+//------------------------------------------------------------------------------------------------------------------------
//
void L1BaseStationInfo::Reset()
{
@@ -133,11 +134,7 @@ void L1BaseStationInfo::Reset()
this->Swap(other);
}
-
-//
-// GETTERS
-//
-//----------------------------------------------------------------------------------------------------------------------//
+//------------------------------------------------------------------------------------------------------------------------
//
const L1Station& L1BaseStationInfo::GetL1Station() const
{
@@ -148,11 +145,27 @@ const L1Station& L1BaseStationInfo::GetL1Station() const
return fL1Station;
}
+//------------------------------------------------------------------------------------------------------------------------
//
-// SETTERS
-//
+const L1Material& L1BaseStationInfo::GetMaterialMap() const
+{
+ if (!fInitController.IsFinalized()) {
+ LOG(fatal) << "L1BaseStationInfo::GetMaterialMap: attempt of getting the material map object from uninitialized "
+ "L1BaseStation class (detectorID = "
+ << static_cast<int>(fDetectorID) << ", stationID = " << fStationID << ")";
+ }
+
+ if (fManagementFlags[static_cast<int>(EManagementFlag::kThicknessMapMoved)]) {
+ LOG(fatal) << "L1BaseStationInfo::GetMaterialMap: attempt of getting the material map, which has been moved. The "
+ "thickness map instance have been "
+ << "already took from the L1BaseStationInfo object (detectorID = " << static_cast<int>(fDetectorID)
+ << ", stationID = " << fStationID << ")";
+ }
-//----------------------------------------------------------------------------------------------------------------------//
+ return fThicknessMap;
+}
+
+//------------------------------------------------------------------------------------------------------------------------
//
void L1BaseStationInfo::SetDetectorID(L1DetectorID inID)
{
@@ -165,7 +178,7 @@ void L1BaseStationInfo::SetDetectorID(L1DetectorID inID)
}
}
-//----------------------------------------------------------------------------------------------------------------------//
+//------------------------------------------------------------------------------------------------------------------------
//
void L1BaseStationInfo::SetRmax(double inRmax)
{
@@ -173,7 +186,7 @@ void L1BaseStationInfo::SetRmax(double inRmax)
fInitController.SetFlag(EInitKey::kRmax);
}
-//----------------------------------------------------------------------------------------------------------------------//
+//------------------------------------------------------------------------------------------------------------------------
//
void L1BaseStationInfo::SetRmin(double inRmin)
{
@@ -181,28 +194,10 @@ void L1BaseStationInfo::SetRmin(double inRmin)
fInitController.SetFlag(EInitKey::kRmin);
}
-//----------------------------------------------------------------------------------------------------------------------//
+//------------------------------------------------------------------------------------------------------------------------
//
-void L1BaseStationInfo::SetFieldSlice(const double* Cx, const double* Cy, const double* Cz)
-{
- if (fInitController.GetFlag(EInitKey::kFieldSlice)) {
- LOG(warn) << "L1BaseStationInfo::SetFieldSlice: Attempt to redifine field slice for station with detectorID = "
- << static_cast<int>(fDetectorID) << " and stationID = " << fStationID << ". Redifinition ignored";
- return;
- }
-
- for (int idx = 0; idx < L1Parameters::kMaxNFieldApproxCoefficients; ++idx) {
- fL1Station.fieldSlice.cx[idx] = Cx[idx];
- fL1Station.fieldSlice.cy[idx] = Cy[idx];
- fL1Station.fieldSlice.cz[idx] = Cz[idx];
- }
-
- fInitController.SetFlag(EInitKey::kFieldSlice);
-}
-
-//----------------------------------------------------------------------------------------------------------------------//
-//
-void L1BaseStationInfo::SetFieldSlice(const std::function<void(const double (&xyz)[3], double (&B)[3])>& getFieldValue)
+void L1BaseStationInfo::SetFieldFunction(
+ const std::function<void(const double (&xyz)[3], double (&B)[3])>& getFieldValue)
{
if (fInitController.GetFlag(EInitKey::kFieldSlice)) {
LOG(warn) << "L1BaseStationInfo::SetFieldSlice: Attempt to redifine field slice for station with detectorID = "
@@ -287,7 +282,7 @@ void L1BaseStationInfo::SetFieldSlice(const std::function<void(const double (&xy
fInitController.SetFlag(EInitKey::kFieldSlice);
}
-//----------------------------------------------------------------------------------------------------------------------//
+//------------------------------------------------------------------------------------------------------------------------
//
void L1BaseStationInfo::SetFieldStatus(int fieldStatus)
{
@@ -295,11 +290,11 @@ void L1BaseStationInfo::SetFieldStatus(int fieldStatus)
fInitController.SetFlag(EInitKey::kFieldStatus);
}
-//----------------------------------------------------------------------------------------------------------------------//
+//------------------------------------------------------------------------------------------------------------------------
//
void L1BaseStationInfo::SetFrontBackStripsGeometry(double frontPhi, double frontSigma, double backPhi, double backSigma)
{
- //----- Original code from L1Algo ---------------------------------------------------------------------//
+ //----- Original code from L1Algo -----------------------------------------------------------------------
double cFront = cos(frontPhi);
double sFront = sin(frontPhi);
double cBack = cos(backPhi);
@@ -327,7 +322,7 @@ void L1BaseStationInfo::SetFrontBackStripsGeometry(double frontPhi, double front
fL1Station.yInfo.cos_phi = cBack / (cBack * sFront - cFront * sBack);
fL1Station.yInfo.sin_phi = -cFront / (cBack * sFront - cFront * sBack);
fL1Station.yInfo.sigma2 = fL1Station.XYInfo.C11;
- //-----------------------------------------------------------------------------------------------------//
+ //-------------------------------------------------------------------------------------------------------
fInitController.SetFlag(EInitKey::kStripsFrontPhi);
fInitController.SetFlag(EInitKey::kStripsFrontSigma);
@@ -335,9 +330,9 @@ void L1BaseStationInfo::SetFrontBackStripsGeometry(double frontPhi, double front
fInitController.SetFlag(EInitKey::kStripsBackSigma);
}
-//----------------------------------------------------------------------------------------------------------------------//
+//------------------------------------------------------------------------------------------------------------------------
//
-void L1BaseStationInfo::SetMaterial(double inThickness, double inRL)
+void L1BaseStationInfo::SetMaterialSimple(double inThickness, double inRL)
{
//L1MASSERT(0, inRL, "Attempt of entering zero inRL (radiational length) value");
@@ -345,11 +340,90 @@ void L1BaseStationInfo::SetMaterial(double inThickness, double inRL)
fL1Station.materialInfo.RL = inRL;
fL1Station.materialInfo.RadThick = fL1Station.materialInfo.thick / fL1Station.materialInfo.RL;
fL1Station.materialInfo.logRadThick = log(fL1Station.materialInfo.RadThick);
- fInitController.SetFlag(EInitKey::kMaterialInfoThick);
- fInitController.SetFlag(EInitKey::kMaterialInfoRL);
+ fInitController.SetFlag(EInitKey::kMaterialInfo);
+}
+
+//------------------------------------------------------------------------------------------------------------------------
+//
+void L1BaseStationInfo::SetMaterialMap(const L1Material& thicknessMap,
+ const std::function<void(L1Material& thicknessMap)>& correction)
+{
+ if (!fInitController.GetFlag(EInitKey::kThicknessMap)) {
+ fThicknessMap = thicknessMap;
+ if (correction) { correction(fThicknessMap); }
+ fInitController.SetFlag(EInitKey::kThicknessMap);
+ }
+ else {
+ LOG(warn) << "L1BaseStationInfo::SetMaterialMap: attempt to reinitialize the material map";
+ }
}
-//----------------------------------------------------------------------------------------------------------------------//
+//------------------------------------------------------------------------------------------------------------------------
+//
+void L1BaseStationInfo::SetMaterialMap(
+ const L1Material& thicknessMap,
+ const std::function<void(L1Material& thicknessMap, const L1MaterialInfo& homogenious)>& correction)
+{
+ if (!fInitController.GetFlag(EInitKey::kMaterialInfo)) {
+ LOG(fatal)
+ << "L1BaseStationInfo::SetMaterialMap: It is impossible to setup the station material thickness map with a "
+ << "correction, which utilizes information on the average station thickness and radiation length. Please, insure "
+ "to "
+ << "set material info with the L1BaseStationInfo::SetMaterialSimple method before setting the thickness map or "
+ << "use correction without the information on average station thickness and radiation length.";
+ }
+
+ if (!fInitController.GetFlag(EInitKey::kThicknessMap)) {
+ fThicknessMap = thicknessMap;
+ correction(fThicknessMap, fL1Station.materialInfo);
+ fInitController.SetFlag(EInitKey::kThicknessMap);
+ }
+ else {
+ LOG(warn) << "L1BaseStationInfo::SetMaterialMap: attempt to reinitialize the material map";
+ }
+}
+
+//------------------------------------------------------------------------------------------------------------------------
+//
+void L1BaseStationInfo::SetMaterialMap(L1Material&& thicknessMap,
+ const std::function<void(L1Material& thicknessMap)>& correction) noexcept
+{
+ if (!fInitController.GetFlag(EInitKey::kThicknessMap)) {
+ fThicknessMap = std::move(thicknessMap);
+ if (correction) { correction(fThicknessMap); }
+ fInitController.SetFlag(EInitKey::kThicknessMap);
+ }
+ else {
+ LOG(warn) << "L1BaseStationInfo::SetMaterialMap: attempt to reinitialize the material map";
+ }
+}
+
+//------------------------------------------------------------------------------------------------------------------------
+//
+void L1BaseStationInfo::SetMaterialMap(
+ L1Material&& thicknessMap,
+ const std::function<void(L1Material& thicknessMap, const L1MaterialInfo& homogenious)>& correction) noexcept
+{
+ if (!fInitController.GetFlag(EInitKey::kMaterialInfo)) {
+ LOG(fatal)
+ << "L1BaseStationInfo::SetMaterialMap: It is impossible to setup the station material thickness map with a "
+ << "correction, which utilizes information on the average station thickness and radiation length. Please, insure "
+ "to "
+ << "set material info with the L1BaseStationInfo::SetMaterialSimple method before setting the thickness map or "
+ << "use correction without the information on average station thickness and radiation length.";
+ }
+
+ if (!fInitController.GetFlag(EInitKey::kThicknessMap)) {
+ fThicknessMap = std::move(thicknessMap);
+ correction(fThicknessMap, fL1Station.materialInfo);
+ fInitController.SetFlag(EInitKey::kThicknessMap);
+ }
+ else {
+ LOG(warn) << "L1BaseStationInfo::SetMaterialMap: attempt to reinitialize the material map";
+ }
+}
+
+//------------------------------------------------------------------------------------------------------------------------
//
void L1BaseStationInfo::SetStationID(int inID)
{
@@ -362,7 +436,7 @@ void L1BaseStationInfo::SetStationID(int inID)
}
}
-//----------------------------------------------------------------------------------------------------------------------//
+//------------------------------------------------------------------------------------------------------------------------
//
void L1BaseStationInfo::SetStationType(int inType)
{
@@ -375,7 +449,7 @@ void L1BaseStationInfo::SetStationType(int inType)
}
}
-//----------------------------------------------------------------------------------------------------------------------//
+//------------------------------------------------------------------------------------------------------------------------
//
void L1BaseStationInfo::SetXmax(double aSize)
{
@@ -383,7 +457,7 @@ void L1BaseStationInfo::SetXmax(double aSize)
fInitController.SetFlag(EInitKey::kXmax);
}
-//----------------------------------------------------------------------------------------------------------------------//
+//------------------------------------------------------------------------------------------------------------------------
//
void L1BaseStationInfo::SetYmax(double aSize)
{
@@ -391,7 +465,7 @@ void L1BaseStationInfo::SetYmax(double aSize)
fInitController.SetFlag(EInitKey::kYmax);
}
-//----------------------------------------------------------------------------------------------------------------------//
+//------------------------------------------------------------------------------------------------------------------------
//
void L1BaseStationInfo::SetTimeInfo(int inTimeInfo)
{
@@ -399,7 +473,7 @@ void L1BaseStationInfo::SetTimeInfo(int inTimeInfo)
fInitController.SetFlag(EInitKey::kTimeInfo);
}
-//----------------------------------------------------------------------------------------------------------------------//
+//------------------------------------------------------------------------------------------------------------------------
//
void L1BaseStationInfo::SetTimeResolution(double dt)
{
@@ -407,7 +481,7 @@ void L1BaseStationInfo::SetTimeResolution(double dt)
fInitController.SetFlag(EInitKey::kTimeResolution);
}
-//----------------------------------------------------------------------------------------------------------------------//
+//------------------------------------------------------------------------------------------------------------------------
//
void L1BaseStationInfo::SetTrackingStatus(bool flag)
{
@@ -415,7 +489,7 @@ void L1BaseStationInfo::SetTrackingStatus(bool flag)
fInitController.SetFlag(EInitKey::kTrackingStatus);
}
-//----------------------------------------------------------------------------------------------------------------------//
+//------------------------------------------------------------------------------------------------------------------------
//
void L1BaseStationInfo::SetZ(double inZ)
{
@@ -424,7 +498,7 @@ void L1BaseStationInfo::SetZ(double inZ)
fInitController.SetFlag(EInitKey::kZ);
}
-//----------------------------------------------------------------------------------------------------------------------//
+//------------------------------------------------------------------------------------------------------------------------
//
void L1BaseStationInfo::Swap(L1BaseStationInfo& other) noexcept
{
@@ -435,10 +509,34 @@ void L1BaseStationInfo::Swap(L1BaseStationInfo& other) noexcept
std::swap(fYmax, other.fYmax);
std::swap(fZPos, other.fZPos);
std::swap(fL1Station, other.fL1Station);
+ std::swap(fThicknessMap, other.fThicknessMap);
std::swap(fInitController, other.fInitController);
}
-//----------------------------------------------------------------------------------------------------------------------//
+//------------------------------------------------------------------------------------------------------------------------
+//
+L1Material&& L1BaseStationInfo::TakeMaterialMap()
+{
+ if (!fInitController.IsFinalized()) {
+ LOG(fatal) << "L1BaseStationInfo::GetMaterialMap: attempt of getting the material map object from uninitialized "
+ "L1BaseStation class (detectorID = "
+ << static_cast<int>(fDetectorID) << ", stationID = " << fStationID << ")";
+ }
+
+ if (fManagementFlags[static_cast<int>(EManagementFlag::kThicknessMapMoved)]) {
+ LOG(fatal) << "L1BaseStationInfo::GetMaterialMap: attempt of taking the material map, which has been moved. The "
+ "thickness map instance have been "
+ << "already took from the L1BaseStationInfo object (detectorID = " << static_cast<int>(fDetectorID)
+ << ", stationID = " << fStationID << ")";
+ }
+ else {
+ fManagementFlags[static_cast<int>(EManagementFlag::kThicknessMapMoved)] = true;
+ }
+
+ return std::move(fThicknessMap);
+}
+
+//------------------------------------------------------------------------------------------------------------------------
//
std::string L1BaseStationInfo::ToString(int verbosityLevel, int indentLevel) const
{
diff --git a/reco/L1/L1Algo/L1BaseStationInfo.h b/reco/L1/L1Algo/L1BaseStationInfo.h
index a080a0476c5b2e45e35468bdcc883cfa666a1a79..21d631f5f21daea966e24c35aa11ed95bbf9d08f 100644
--- a/reco/L1/L1Algo/L1BaseStationInfo.h
+++ b/reco/L1/L1Algo/L1BaseStationInfo.h
@@ -19,6 +19,7 @@
#include "L1ObjectInitController.h"
#include "L1Station.h"
// C++ std
+#include <bitset>
#include <functional>
#include <string>
@@ -26,6 +27,14 @@ enum class L1DetectorID;
/// A base class which provides interface to L1Algo station geometry
class L1BaseStationInfo {
+ /// Enumeration for internal logic control
+ enum class EManagementFlag
+ {
+ kThicknessMapMoved, ///< if the thickness map was moved from the L1BaseStationInfo instance
+ kEnd
+ };
+ using ManagementFlags_t = std::bitset<static_cast<int>(EManagementFlag::kEnd)>;
+
public:
/// Enumeration of fields, which must be initialized so the object can pass the threshold
enum class EInitKey
@@ -43,8 +52,8 @@ public:
kZ, ///< z coordinate of the station position
kRmin, ///< internal radius of station (gap size)
kRmax, ///< exteranl radius of station
- kMaterialInfoThick, ///< thickness of the station
- kMaterialInfoRL, ///< rad length of the station
+ kMaterialInfo, ///< thickness and rad. length of the station
+ kThicknessMap, ///< thickness map of the station (optional?)
kFieldSlice, ///< L1Station.L1FieldSlice object initialization
kStripsFrontPhi, ///< strips geometry initialization
kStripsFrontSigma, ///<
@@ -55,86 +64,112 @@ public:
kEnd
};
- using L1ObjectInitController_t = L1ObjectInitController<static_cast<int>(EInitKey::kEnd), EInitKey>;
+ using InitController_t = L1ObjectInitController<static_cast<int>(EInitKey::kEnd), EInitKey>;
//
// CONSTRUCTORS AND DESTRUCTORS
//
/// Default constructor
L1BaseStationInfo() noexcept;
+
/// Constructor from stationID and detectorID
L1BaseStationInfo(L1DetectorID detetorID, int stationID) noexcept;
+
/// Destructor
~L1BaseStationInfo() noexcept;
+
/// Copy constructor
L1BaseStationInfo(const L1BaseStationInfo& other) noexcept;
+
/// Move constructor
L1BaseStationInfo(L1BaseStationInfo&& other) noexcept;
+
/// Copy assignment operator
L1BaseStationInfo& operator=(const L1BaseStationInfo& other) noexcept;
+
/// Move assignment operator
L1BaseStationInfo& operator=(L1BaseStationInfo&& other) noexcept;
- //
- // BASIC METHODS
- //
- /// Less operator for L1BaseStationInfo object to perform their sorts. Sorting is carried out first by fDetectorID,
- /// and then by fStationID
+ /// Less operator for L1BaseStationInfo object to perform sorts. Sorting is carried out first by fDetectorID,
+ /// and then by fStationID. The detectorID must be defined in L1DetectorID enum class
bool operator<(const L1BaseStationInfo& right) const
{
return fDetectorID != right.fDetectorID ? fDetectorID < right.fDetectorID : fStationID < right.fStationID;
}
- //
- // GETTERS
- //
/// Gets detector ID
L1DetectorID GetDetectorID() const { return fDetectorID; }
+
/// Gets a coefficient with idx for the field x-component approximation
fvec GetFieldSliceCx(int idx) const { return fL1Station.fieldSlice.cx[idx]; }
+
/// Gets a coefficient with idx for the field y-component approximation
fvec GetFieldSliceCy(int idx) const { return fL1Station.fieldSlice.cy[idx]; }
+
/// Gets a coefficient with idx for the field z-component approximation
fvec GetFieldSliceCz(int idx) const { return fL1Station.fieldSlice.cz[idx]; }
+
/// Gets array of the coefficients for the field x-component approximation
const fvec* GetFieldSliceCx() const { return fL1Station.fieldSlice.cx; }
+
/// Gets array of the coefficients for the field y-component approximation
const fvec* GetFieldSliceCy() const { return fL1Station.fieldSlice.cy; }
+
/// Gets array of the coefficients for the field z-component approximation
const fvec* GetFieldSliceCz() const { return fL1Station.fieldSlice.cz; }
+
/// Gets field status: 0 - station is outside the field, 1 - station is inside the field
int GetFieldStatus() const { return fL1Station.fieldStatus; }
+
/// Gets a const reference to the L1ObjectInitController object
- const L1ObjectInitController_t& GetInitController() const { return fInitController; }
+ const InitController_t& GetInitController() const { return fInitController; }
+
/// Gets a reference to L1Station info field of the L1BaseStation info
const L1Station& GetL1Station() const;
- /// Gets material thickness
- fvec GetMaterialThickness() const { return fL1Station.materialInfo.thick; }
+
+ /// Gets a reference to L1Material map
+ const L1Material& GetMaterialMap() const;
+
+ /// Gets station thickness
+ fvec GetThickness() const { return fL1Station.materialInfo.thick; }
+
/// Gets the radiation length of the station material
fvec GetMaterialRadLength() const { return fL1Station.materialInfo.RL; }
- /// Gets the relative material thickness in units of the radiational length
- fvec GetMaterialRadThick() const { return fL1Station.materialInfo.RadThick; }
- /// Gets log of the relative material thickness in units of the radiational length
- fvec GetMaterialLogRadThick() const { return fL1Station.materialInfo.logRadThick; }
- /// Gets min transverse size of the station
+
+ /// Gets the station thickness in units of the radiation length
+ fvec GetRadThick() const { return fL1Station.materialInfo.RadThick; }
+
+ /// Gets log of the station thickness in units of the radiation length
+ fvec GetLogRadThick() const { return fL1Station.materialInfo.logRadThick; }
+
+ /// Gets min transverse size of the station [cm]
fvec GetRmin() const { return fL1Station.Rmin; }
- /// Gets max transverse size of the station
+
+ /// Gets max transverse size of the station [cm]
fvec GetRmax() const { return fL1Station.Rmax; }
+
/// Gets station ID
int GetStationID() const { return fStationID; }
+
/// Gets station type
int GetStationType() const { return fL1Station.type; }
+
/// Gets time resolution
fvec GetTimeResolution() const { return fL1Station.dt; }
+
/// Gets tracking status: true - station is active for tracking, false - station exists, but not used in tracking
bool GetTrackingStatus() const { return fTrackingStatus; }
+
/// Gets maximum distance between station center and its edge in x direction
double GetXmax() const { return fXmax; }
+
/// Gets maximum distance between station center and its edge in y direction
double GetYmax() const { return fYmax; }
- /// Gets double precised z position of the station
+
+ /// Gets double precised z position of the station [cm]
double GetZdouble() const { return fZPos; }
- /// Gets SIMD vectorized z position of the station
+
+ /// Gets SIMD vectorized z position of the station [cm]
fvec GetZsimdVec() const { return fL1Station.z; }
/// Prints registered fields
@@ -142,73 +177,115 @@ public:
/// verbosity = 1: print basic L1Station fields
/// verbosity = 2: print all L1Staion fields
void Print(int verbosity = 0) const;
+
/// Resets fields to the default values
void Reset();
- //
- // SETTERS
- //
/// Sets detector ID
void SetDetectorID(L1DetectorID inID);
+
/// Sets flag: true - station is placed in field, false - station is placed outside the field
void SetFieldStatus(int fieldStatus);
- /// Sets field approximation at the station plain
- /// \param Cx Array of approximation coefficients for x field component
- /// \param Cy Array of approximation coefficients for y field component
- /// \param Cz Array of approximation coefficients for z field component
- void SetFieldSlice(const double* Cx, const double* Cy, const double* Cz);
+
/// Sets arrays of the approcimation
/// \param getField A user function, which gets a xyz array of position coordinates and fills B array
/// of magnetic field components in position
- void SetFieldSlice(const std::function<void(const double (&xyz)[3], double (&B)[3])>& getFieldValue);
+ void SetFieldFunction(const std::function<void(const double (&xyz)[3], double (&B)[3])>& getFieldValue);
+
/// Sets stereo angles and sigmas for front and back strips
- /// \param f_phi Stereoangle of front strips
- /// \param f_sigma Sigma of front strips
- /// \param b_phi Stereoangle of back strips
- /// \param b_sigma Sigma of back strips
+ /// \param f_phi Stereoangle of front strips [rad]
+ /// \param f_sigma Sigma of front strips [rad]
+ /// \param b_phi Stereoangle of back strips [rad]
+ /// \param b_sigma Sigma of back strips [rad]
void SetFrontBackStripsGeometry(double fPhi, double fSigma, double bPhi, double bSigma);
- /// Sets station thickness and radiational length
- /// \param thickness thickness of station
- /// \param radiationalLength radiational length of station
- void SetMaterial(double thickness, double radiationalLength);
- /// Sets max transverse size of the station
+
+ /// Sets station thickness and radiation length
+ /// \param thickness Thickness of station [arb. units]
+ /// \param radiationLength Radiation length of station [arb. units]
+ void SetMaterialSimple(double thickness, double radiationLength);
+
+ /// Sets station thickness in units of radiation length mapped vs. position in XY plane (copy semantics)
+ /// \param thicknessMap Map of station thickness in units of radiation length
+ /// \param correction User corrections to material map: thicknessMap - reference to the material map
+ void SetMaterialMap(const L1Material& thicknessMap,
+ const std::function<void(L1Material& thicknessMap)>& correction = nullptr);
+
+ /// Sets station thickness in units of radiation length mapped vs. position in XY plane (copy semantics)
+ /// \param thicknessMap Map of station thickness in units of radiation length
+ /// \param correction User corrections to material map: thicknessMap - reference to the material map,
+ /// homogenious - object, which keeps expected homogenious values of station thickness and rad. length
+ void
+ SetMaterialMap(const L1Material& thicknessMap,
+ const std::function<void(L1Material& thicknessMap, const L1MaterialInfo& homogenious)>& correction);
+
+ /// Sets station thickness in units of radiation length mapped vs. position in XY plane (move semantics)
+ /// \param thicknessMap Map of station thickness in units of radiation length
+ /// \param correction User corrections to material map: thicknessMap - reference to the material map
+ void SetMaterialMap(L1Material&& thicknessMap,
+ const std::function<void(L1Material& thicknessMap)>& correction = nullptr) noexcept;
+
+ /// Sets station thickness in units of radiation length mapped vs. position in XY plane (move semantics)
+ /// \param thicknessMap Map of station thickness in units of radiation length
+ /// \param correction User corrections to material map: thicknessMap - reference to the material map,
+ /// homogenious - object, which keeps expected homogenious values of station thickness and rad. length
+ void SetMaterialMap(
+ L1Material&& thicknessMap,
+ const std::function<void(L1Material& thicknessMap, const L1MaterialInfo& homogenious)>& correction) noexcept;
+
+ /// Sets max transverse size of the station []
void SetRmax(double inRmax);
- /// Sets min transverse size of the station
+
+ /// Sets min transverse size of the station []
void SetRmin(double inRmin);
+
/// Sets station ID
void SetStationID(int inID);
+
/// Sets type of station
void SetStationType(int inType); // TODO: this is a temporary solution (S.Zh.)
+
/// Sets flag: 0 - time information is not provided by this detector type
/// 1 - time information is provided by the detector and can be used in tracking
void SetTimeInfo(int inTimeInfo);
+
/// Sets time resolution
void SetTimeResolution(double dt);
+
/// Sets tracking status: true - station is active for tracking, false - station exists, but not used in tracking
void SetTrackingStatus(bool flag);
+
/// Sets maximum distance between station center and its edge in x direction
void SetXmax(double aSize);
+
/// Sets maximum distance between station center and its edge in y direction
void SetYmax(double aSize);
+
/// Sets nominal z position of the station
void SetZ(double inZ);
- /// Swap method for easier implementation of above ones (NOTE: all the fields must be accounted carefully here)
+ /// Swap method for easy implementation of move constructor and copy and move assignment operator
void Swap(L1BaseStationInfo& other) noexcept;
+ /// Takes an instance of L1Material map
+ /// This method can be called only once. It will catch the content of L1Material map (using move semantics), and after the warning will be generated
+ /// both for GetMaterialMap and TakeMaterialMap
+ L1Material&& TakeMaterialMap();
+
/// String representation of class contents
- /// \param indentLevel number of indent characters in the output
+ /// \param indentLevel Number of indent characters in the output
std::string ToString(int verbosityLevel = 0, int indentLevel = 0) const;
private:
L1DetectorID fDetectorID {static_cast<L1DetectorID>(0)}; ///< Detector ID
int fStationID {-1}; ///< Station ID
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
- double fZPos {0}; ///< z position of the station in double precision, used in field approximation
- L1Station fL1Station {}; ///< L1Station structure, describes a station in L1Algo
- L1ObjectInitController_t fInitController {}; ///< Class fileds initialization flags
+ 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
+ double fZPos {0}; ///< z position of the station in double precision, used in field approximation
+ L1Station fL1Station {}; ///< L1Station structure, describes a station in L1Algo
+ L1Material fThicknessMap {}; ///< Map of station thickness in units of radiation length
+ InitController_t fInitController {}; ///< Class fileds initialization flags
+ ManagementFlags_t fManagementFlags {}; ///< bitset flags to manage internal behaviour of the class
};
/// swap function for two L1BaseStationInfo objects, expected to be used instead of std::swap
diff --git a/reco/L1/L1Algo/L1InitManager.cxx b/reco/L1/L1Algo/L1InitManager.cxx
index 6c5335391a4aa9ae3843f483a4ce6889983db3f2..be625c35462b23f2a08c79e0a30344249096e14a 100644
--- a/reco/L1/L1Algo/L1InitManager.cxx
+++ b/reco/L1/L1Algo/L1InitManager.cxx
@@ -41,7 +41,19 @@ void L1InitManager::AddStation(const L1BaseStationInfo& inStation)
if (isStationActive) {
// initialize magnetic field slice
L1BaseStationInfo inStationCopy = L1BaseStationInfo(inStation); // make a copy of station so it can be initialized
- inStationCopy.SetFieldSlice(fFieldFunction);
+ inStationCopy.SetFieldFunction(fFieldFunction);
+
+ // check, if material map is used
+ if (!inStationCopy.GetInitController().GetFlag(L1BaseStationInfo::EInitKey::kThicknessMap)) {
+ LOG(warn) << "Station material map was not set for detectorID = "
+ << static_cast<int>(inStationCopy.GetDetectorID()) << ", stationID = " << inStationCopy.GetStationID()
+ << ". Homogenious material budget will be used: " << inStationCopy.GetRadThick()[0];
+ L1Material material;
+ material.SetBins(1, 100);
+ material.SetRadThick(0, 0, inStationCopy.GetRadThick()[0]);
+ inStationCopy.SetMaterialMap(std::move(material));
+ }
+
// Check station init
{
LOG(debug) << "L1InitManager::AddStation:(original) L1BaseStationInfo "
@@ -313,12 +325,47 @@ void L1InitManager::TransferL1StationArray(std::array<L1Station, L1Parameters::k
auto destinationArrayIterator = destinationArray.begin();
for (const auto& item : fStationsInfo) {
- *destinationArrayIterator = std::move(item.GetL1Station());
+ *destinationArrayIterator = item.GetL1Station();
++destinationArrayIterator;
}
- LOG(info) << "L1InitManager: L1Station vector was successfully transfered to L1Algo core :)";
+ LOG(info) << "L1InitManager: L1Station vector was successfully transfered to L1Algo core";
}
+//-----------------------------------------------------------------------------------------------------------------------
+//
+void L1InitManager::TransferL1MaterialArray(L1Vector<L1Material>& destinationArray)
+{
+ //
+ // 1) Check, if all fields of this were initialized
+ //
+ {
+ std::stringstream aStream;
+ aStream << "Attempt to pass L1Station array to L1Algo core before all necessary fields initialization\n"
+ << "L1InitManager " << fInitController.ToString();
+ L1MASSERT(0, fInitController.IsFinalized(), aStream.str().c_str());
+ }
+ //
+ // 2) Check, if destinationArraySize is enough for the transfer
+ //
+ {
+ int nStationsTotal = this->GetNstationsActive();
+ std::stringstream aStream;
+ aStream << "Destination array size (" << destinationArray.capacity()
+ << ") is smaller then the actual number of active tracking stations (" << nStationsTotal << ")";
+ L1MASSERT(0, nStationsTotal <= static_cast<int>(destinationArray.capacity()), aStream.str().c_str());
+ }
+
+ auto destinationArrayIterator = destinationArray.begin();
+ for (auto it = fStationsInfo.begin(); it != fStationsInfo.end(); ++it) {
+ auto node = fStationsInfo.extract(it);
+ *destinationArrayIterator = std::move(node.value().TakeMaterialMap());
+ fStationsInfo.insert(std::move(node));
+ ++destinationArrayIterator;
+ }
+ LOG(info) << "L1InitManager: L1Material vector was successfully transfered to L1Algo core";
+}
+
+
//
// INIT CHECKERS
//
diff --git a/reco/L1/L1Algo/L1InitManager.h b/reco/L1/L1Algo/L1InitManager.h
index 456fcfc40f7cfab814ab4d12efb867d4df847b54..47191c1a6407af06d713454ed40bc9442c766dcb 100644
--- a/reco/L1/L1Algo/L1InitManager.h
+++ b/reco/L1/L1Algo/L1InitManager.h
@@ -89,34 +89,37 @@ private:
using L1DetectorIDIntMap_t = std::unordered_map<L1DetectorID, int, L1Utils::EnumClassHash>;
using L1DetectorIDSet_t = std::set<L1DetectorID>;
using L1FieldFunction_t = std::function<void(const double (&xyz)[3], double (&B)[3])>;
- using L1ObjectInitController_t = L1ObjectInitController<static_cast<int>(EInitKey::kEnd), EInitKey>;
+ using InitController_t = L1ObjectInitController<static_cast<int>(EInitKey::kEnd), EInitKey>;
public:
- //
- // CONSTRUCTORS AND DESTRUCTOR
- //
/// Default constructor
L1InitManager() = delete;
+
/// Constructor from ptr to L1Paramters object
L1InitManager(L1Parameters* pParameters);
+
/// Destructor
~L1InitManager() = default;
+
/// Copy constructor is forbidden
L1InitManager(const L1InitManager& /*other*/) = delete;
+
/// Move constructor is forbidden
L1InitManager(L1InitManager&& /*other*/) = delete;
+
/// Copy assignment operator is forbidden
L1InitManager& operator=(const L1InitManager& /*other*/) = delete;
+
/// Move assignment operator is forbidden
L1InitManager& operator=(L1InitManager&& /*other*/) = delete;
- //
- // BASIC METHODS
- //
+
/// Adds another station of a given type using reference to a L1BaseStationInfo object
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); }
@@ -126,33 +129,36 @@ public:
// NOTE: This method calls checkers of large fields initializations like a station or an iteration. The method must be
// called in the L1Algo class. (S.Zharko)
- //
- // GETTERS
- //
- /// Gets a set of active detectors for this analysis
- const L1DetectorIDSet_t& GetActiveDetectorIDs() const { return fActiveDetectorIDs; }
/// Gets ghost suppression flag
int GetGhostSuppression() const { return fGhostSuppression; }
+
/// Gets momentum cutoff
float GetMomentumCutOff() const { return fMomentumCutOff; }
+
/// Gets a const reference to L1ObjectInitController
- const L1ObjectInitController_t& GetInitController() const { return fInitController; }
+ const InitController_t& GetInitController() const { return fInitController; }
+
/// Gets a pointer to L1Parameters instance with a posibility of its fields modification
const L1Parameters* GetParameters() const { return fpParameters; }
+
/// Gets total number of active stations
int GetNstationsActive() const { return fNstationsActive[fNstationsActive.size() - 1]; }
+
/// Gets number of active stations for given detector ID
int GetNstationsActive(L1DetectorID detectorID) const
{
return fNstationsActive[static_cast<L1DetectorID_t>(detectorID)];
}
+
/// Gets total number of stations, provided by setup geometry
int GetNstationsGeometry() const { return fNstationsGeometry[fNstationsGeometry.size() - 1]; }
+
/// Gets number of stations, provided by setup geometry for given detector ID
int GetNstationsGeometry(L1DetectorID detectorID) const
{
return fNstationsGeometry[static_cast<L1DetectorID_t>(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
@@ -162,6 +168,7 @@ public:
+ std::accumulate(fNstationsGeometry.cbegin(), fNstationsGeometry.cbegin() + static_cast<int>(detectorID),
0);
}
+
/// 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
@@ -172,10 +179,13 @@ public:
/// Gets L1FieldRegion object at primary vertex
const L1FieldRegion& GetTargetFieldRegion() const { return fTargetFieldRegion; }
+
/// Gets L1FieldValue object at primary vertex
const L1FieldValue& GetTargetFieldValue() const { return fTargetFieldValue; }
+
/// Gets a target position
const std::array<double, 3>& GetTargetPosition() const { return fTargetPos; }
+
/// Gets tracking level
int GetTrackingLevel() const { return fTrackingLevel; }
@@ -186,56 +196,67 @@ public:
/// Prints a list of CA track finder iterations
void PrintCAIterations(int verbosityLevel = 0) const;
+
/// Prints a list of stations
void PrintStations(int verbosityLevel = 0) const;
/// Pushes an CA track finder iteration into a sequence of iteration using reference
void PushBackCAIteration(const L1CAIteration& iteration);
+
/// Pushes an CA track finder iteration into a sequence of iteration using raw pointer
void PushBackCAIteration(const L1CAIteration* pIteration) { PushBackCAIteration(*pIteration); }
+
/// Pushes an CA track finder iteration into a sequence of iteration using std::unique_ptr
void PushBackCAIteration(const std::unique_ptr<L1CAIteration>& puIteration) { PushBackCAIteration(*puIteration); }
- //
- // SETTERS
- //
/// Sets a set of active tracking detector IDs
void SetActiveDetectorIDs(const L1DetectorIDSet_t& detectorIDs);
+
/// Sets a number of CA track finder iterations to provide initialization cross-check
void SetCAIterationsNumberCrosscheck(int nIterations);
+
/// Sets a magnetic field function, which will be applied for all the stations
void SetFieldFunction(const L1FieldFunction_t& fieldFcn);
+
///
void SetGhostSuppression(int ghostSuppression);
+
///
void SetMomentumCutOff(float momentumCutOff);
+
///
void SetTrackingLevel(int trackingLevel);
+
/// Sets a number of actual stations for a particular tracking detector ID to provide initialization cross-check
void SetNstations(L1DetectorID detectorID, int nStations);
+
/// Sets target poisition
void SetTargetPosition(double x, double y, double z);
/// Transfers an array of L1Stations formed inside a set of L1BaseStationInfo to a destination std::array
void TransferL1StationArray(std::array<L1Station, L1Parameters::kMaxNstations>& destinationArray);
+ /// Transfers an array of L1Material tables formed inside a set of L1BaseStationInfo to a destination std::array
+ void TransferL1MaterialArray(L1Vector<L1Material>& destinationVector);
+
private:
/// Checker for L1CAIteration container initialization (sets EInitKey::kCAIterations)
/// \return true If all L1CAIteration objects were initialized properly
void CheckCAIterationsInit();
+
/// Checker for L1BaseStationInfo set initialization (sets EInitKey::kStationsInfo)
/// \return true If all L1BaseStationInfo objects were initialized properly. Similar effect can be achieved by
void CheckStationsInfoInit();
/* Basic fields */
- L1ObjectInitController_t fInitController {}; ///< Initialization flags
+ InitController_t fInitController {}; ///< Initialization flags
L1DetectorIDSet_t fActiveDetectorIDs {}; ///< Set of tracking detectors, active during this analysis session
/* Target fields */
- std::array<double, /*nDimensions=*/3> fTargetPos {}; ///< Nominal target position coordinates
+ std::array<double, /*nDimensions=*/3> fTargetPos {}; ///< Nominal target position coordinates [cm]
/* Stations related fields */
@@ -244,9 +265,11 @@ private:
/// 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 subystem. The last array element corresponds to the total number of stations.
std::array<int, L1Parameters::kMaxNdetectors + 1> fNstationsActive {};
+
/// Actual numbers of stations, provided by geometry. Index of an array element (except the last one) corresponds to a given
/// L1DetectorID of the detector subystem. The last array element corresponds to the total number of stations.
std::array<int, L1Parameters::kMaxNdetectors + 1> fNstationsGeometry {};
+
/// Map of the actual detector indeces to the active detector indeces
/// 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:
@@ -256,8 +279,7 @@ private:
/// A function which returns magnetic field vector B in a radius-vector xyz
L1FieldFunction_t fFieldFunction {[](const double (&)[3], double (&)[3]) {}};
- // NOTE: Stations of daetectors which will not be assigned as active, will not be included in the tracking!!!!!!!
- // NOTE: fTotalNumberOfStations is excess field for logic, but it's imortant to track L1Algo initialization
+ // NOTE: Stations of daetectors which will not be assigned as active, will not be included in the tracking!
/* Vertex related fields */
@@ -266,12 +288,11 @@ private:
/* CA track finder iterations related */
- //L1Vector<L1CAIteration> fCAIterationsContainer {}; ///< Container for CA track finder iterations
int fCAIterationsNumberCrosscheck {-1}; ///< Number of iterations to be passed (must be used for cross-checks)
- /// Pointer to L1Parameters object
- // NOTE: Owner of the object is L1Algo instance
- L1Parameters* fpParameters {nullptr};
+
+ L1Parameters* fpParameters {
+ nullptr}; ///< Pointer to L1Parameters object. NOTE: Owner of the object is L1Algo instance
int fTrackingLevel {0}; ///< tracking level
int fGhostSuppression {0}; ///< flag: if true, ghost tracks are suppressed
diff --git a/reco/L1/L1Algo/L1MaterialInfo.cxx b/reco/L1/L1Algo/L1MaterialInfo.cxx
index 6c778f20fab3cd8f8dc586f9c8d78ad2b8efe7ac..3eccb8d896e3f5d01fe6ad444aa4c69495e4eecc 100644
--- a/reco/L1/L1Algo/L1MaterialInfo.cxx
+++ b/reco/L1/L1Algo/L1MaterialInfo.cxx
@@ -8,15 +8,123 @@
#include <sstream>
#include <vector>
+/************************
+ * L1MaterialInfo class *
+ ************************/
+
std::string L1MaterialInfo::ToString(int indentLevel) const
{
std::stringstream aStream {};
// TODO: possibly it is better to place indentChar into L1Parameters (S.Zharko)
constexpr char indentChar = '\t';
std::string indent(indentLevel, indentChar);
- aStream << indent << "Thickness (X) [??]: " << std::setw(12) << std::setfill(' ') << thick[0] << '\n';
- aStream << indent << "Radiational length (X0) [??]:" << std::setw(12) << std::setfill(' ') << RL[0] << '\n';
+ aStream << indent << "Station thickness (X) [??]: " << std::setw(12) << std::setfill(' ') << thick[0] << '\n';
+ aStream << indent << "Radiation length (X0) [??]: " << std::setw(12) << std::setfill(' ') << RL[0] << '\n';
aStream << indent << "X/X0: " << std::setw(12) << std::setfill(' ') << RadThick[0] << '\n';
aStream << indent << "log(X/X0): " << std::setw(12) << std::setfill(' ') << logRadThick[0];
return aStream.str();
}
+
+/********************
+ * L1Material class *
+ ********************/
+
+L1Material::L1Material() {}
+
+//------------------------------------------------------------------------------------------------------------------------------------
+//
+L1Material::~L1Material() noexcept {}
+
+//------------------------------------------------------------------------------------------------------------------------------------
+//
+L1Material::L1Material(const L1Material& other)
+ : fNbins(other.fNbins)
+ , fRmax(other.fRmax)
+ , fFactor(other.fFactor)
+ , fTable(other.fTable)
+{
+ std::cout << "L1Material copy constructor was called\n";
+}
+
+//------------------------------------------------------------------------------------------------------------------------------------
+//
+
+L1Material& L1Material::operator=(const L1Material& other)
+{
+ std::cout << "L1Material copy assignment operator was called\n";
+ if (this != &other) {
+ fNbins = other.fNbins;
+ fRmax = other.fRmax;
+ fFactor = other.fFactor;
+ fTable.clear();
+ fTable.resize(other.fTable.size());
+ for (size_t i = 0; i < other.fTable.size(); ++i) {
+ fTable[i] = other.fTable[i];
+ }
+ }
+ return *this;
+}
+
+//------------------------------------------------------------------------------------------------------------------------------------
+//
+L1Material::L1Material(L1Material&& other) noexcept
+{
+ std::cout << "L1Material move constructor was called\n";
+ this->Swap(other);
+}
+
+//------------------------------------------------------------------------------------------------------------------------------------
+//
+L1Material& L1Material::operator=(L1Material&& other) noexcept
+{
+ std::cout << "L1Material move assignment operator was called\n";
+ if (this != &other) {
+ L1Material tmp(std::move(other));
+ this->Swap(tmp);
+ }
+ return *this;
+}
+
+//------------------------------------------------------------------------------------------------------------------------------------
+//
+float L1Material::GetRadThick(float x, float y) const
+{
+ x = (x < fRmax && x >= -fRmax) ? x : 0;
+ y = (y < fRmax && y >= -fRmax) ? y : 0;
+ int i = static_cast<int>((x + fRmax) * fFactor);
+ int j = static_cast<int>((y + fRmax) * fFactor);
+ i = (i < fNbins && i >= 0) ? i : fNbins / 2;
+ j = (j < fNbins && j >= 0) ? j : fNbins / 2;
+ return fTable[i + j * fNbins];
+}
+
+//------------------------------------------------------------------------------------------------------------------------------------
+//
+fvec L1Material::GetRadThick(fvec x, fvec y) const
+{
+ fvec r;
+ for (int i = 0; i < fvecLen; i++)
+ r[i] = GetRadThick(x[i], y[i]);
+ return r;
+}
+
+//------------------------------------------------------------------------------------------------------------------------------------
+//
+void L1Material::SetBins(int nBins, float stationSize)
+{
+ fNbins = nBins;
+ fRmax = stationSize;
+ fFactor = 0.5 * fNbins / fRmax;
+
+ fTable.resize(fNbins * fNbins);
+}
+
+//------------------------------------------------------------------------------------------------------------------------------------
+//
+void L1Material::Swap(L1Material& other) noexcept
+{
+ std::swap(fNbins, other.fNbins);
+ std::swap(fRmax, other.fRmax);
+ std::swap(fFactor, other.fFactor);
+ std::swap(fTable, other.fTable); // Probably can cause segmentation violation (did not understand)
+}
diff --git a/reco/L1/L1Algo/L1MaterialInfo.h b/reco/L1/L1Algo/L1MaterialInfo.h
index c8e15067826e3375ed49eb55acf453ac470b9cd8..762d048281b59295d49980df115ad98990e734d8 100644
--- a/reco/L1/L1Algo/L1MaterialInfo.h
+++ b/reco/L1/L1Algo/L1MaterialInfo.h
@@ -12,11 +12,13 @@
#include "L1Def.h"
-class L1MaterialInfo {
-public:
- fvec thick {0};
- fvec RL {0};
- fvec RadThick {0};
+/// Class L1MaterialInfo contains SIMDized vector fields of the
+/// The fields of the structure should ONLY be initialized within L1BaseStationInfo::SetMaterial(double, double) method, when the
+/// stations sequence is initialized
+struct L1MaterialInfo {
+ fvec thick {0}; ///< Average thickness of the station in arbitary length units
+ fvec RL {0}; ///< Average radiation length (X0) of the station material in THE SAME UNITS as the thickness
+ fvec RadThick {0}; ///< Average thickness in units of radiation length (X/X0)
fvec logRadThick {0};
/// String representation of class contents
@@ -24,52 +26,67 @@ public:
std::string ToString(int indentLevel = 0) const;
} _fvecalignment;
-
-// TODO: Probably one should refactor this class:
-// 1) Make MBins and RMax private
-// 2) Replace std::vector<std::vector<float>> with std::vector<float> (using recalculation of (i,j) into one idx)
-// (S.Zharko)
+/// Class L1Material describes a map of station thickness in units of radiation length (X0) to the specific point in XY plane
class L1Material {
public:
- L1Material() : table(0), NBins(0), RMax(0.), iD(0.) {};
-
- std::vector<std::vector<float>> table;
-
- // static const int NBins = 100; // TODO file?
- // static const float RMax = 60.f;
- // static const float iD = 0.5*NBins/60.f;//RMax!;
-
- void SetBins(int n, float r)
- {
- NBins = n;
- RMax = r;
- iD = 0.5 * NBins / RMax;
- }
-
- float GetRadThick(float x, float y)
- {
- x = (x < RMax && x >= -RMax) ? x : 0;
- y = (y < RMax && y >= -RMax) ? y : 0;
- int i = static_cast<int>((x + RMax) * iD);
- int j = static_cast<int>((y + RMax) * iD);
- i = (i < NBins && i >= 0) ? i : NBins / 2;
- j = (j < NBins && j >= 0) ? j : NBins / 2;
- return table[i][j];
- }
-
-
- fvec GetRadThick(fvec x, fvec y)
- {
- fvec r;
- for (int i = 0; i < fvecLen; i++)
- r[i] = GetRadThick(x[i], y[i]);
- return r;
- }
-
- int NBins;
- float RMax;
- float iD;
-
+ /// Default constructor
+ L1Material();
+
+ /// Copy constructor
+ L1Material(const L1Material& other);
+
+ /// Copy assignment operator
+ L1Material& operator=(const L1Material& other);
+
+ /// Move constructor
+ L1Material(L1Material&& other) noexcept;
+
+ /// Move assignment operator
+ L1Material& operator=(L1Material&& other) noexcept;
+
+ /// Destructor
+ ~L1Material() noexcept;
+
+ /// Gets number of bins (rows or columns) of the material table
+ int GetNbins() const { return fNbins; }
+
+ /// Gets value of X/X0 in a given cell of the material table by the indeces of the row and column
+ /// \param iBinX Index of table column
+ /// \param iBinY Index of table row
+ float GetRadThick(int iBinX, int iBinY) const { return fTable[iBinX + fNbins * iBinY]; }
+
+ /// Gets material thickness in units of X0 in (x,y) point of the station
+ /// \param x X coordinate of the point [cm]
+ /// \param y Y coordinate of the point [cm]
+ float GetRadThick(float x, float y) const;
+
+ /// Gets material thickness in units of X0 in (x,y) point of the station
+ /// \param x X coordinate of the point [cm] (SIMDized vector)
+ /// \param y Y coordinate of the point [cm] (SIMDized veotor)
+ fvec GetRadThick(fvec x, fvec y) const;
+
+ /// Sets value of material thickness in units of X0 for a given cell of the material table
+ /// WARNING: Indeces of rows and columns in the table runs from 0 to nBins-1 inclusively, where nBins is the number both of rows
+ /// and columns. One should be carefull reading and storing the table from ROOT-file, because iBinX = 0 and iBinY = 0 in the
+ /// SetBinContent method of a ROOT histogram usually defines underflow bin
+ /// \param iBinX Index of table column
+ /// \param iBinY Index of table row
+ /// \param thickness Thickness of the material in units of X0
+ void SetRadThick(int iBinX, int iBinY, float thickness) { fTable[iBinX + fNbins * iBinY] = thickness; }
+
+ /// Sets properties of the material table -- number of rows or columnts and the size of station in XY plane
+ /// \param nBins Number of rows or columns
+ /// \param stationSize Size of station in x and y dimensions [cm]
+ void SetBins(int nBins, float stationSize);
+
+ /// Swap method
+ void Swap(L1Material& other) noexcept;
+
+private:
+ int fNbins {0}; ///< Number of rows (columns) in the material budget table
+ float fRmax {0.f}; ///< Size of the station in x and y dimensions [cm]
+ float fFactor {0.f}; ///< Factor used in the recalculation of point coordinates to row/column id
+ std::vector<float> fTable {}; ///< Material budget table
} _fvecalignment;
#endif
diff --git a/reco/L1/ParticleFinder/CbmL1PFFitter.cxx b/reco/L1/ParticleFinder/CbmL1PFFitter.cxx
index 4db6ce2b8ca225393cc28ce418be30b789249a7f..f8afa2a532c396897acd3dd3cb9a56ce4f7c2128 100644
--- a/reco/L1/ParticleFinder/CbmL1PFFitter.cxx
+++ b/reco/L1/ParticleFinder/CbmL1PFFitter.cxx
@@ -257,13 +257,13 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo)
fvec wIn = (ONE & (initialised));
L1Extrapolate(T, z[i], qp0, fld, &w1);
- if (i == NMvdStations) {
+ if (i == NMvdStations) { // TODO: How a hit can be also a station? (S.Zharko)
fit.L1AddPipeMaterial(T, qp0, wIn);
fit.EnergyLossCorrection(T, fit.PipeRadThick, qp0, fvec(-1.f), wIn);
}
if constexpr (L1Parameters::kIfUseRadLengthTable) {
- fit.L1AddMaterial(T, CbmL1::Instance()->algo->fRadThick[i].GetRadThick(T.x, T.y), qp0, wIn);
- fit.EnergyLossCorrection(T, CbmL1::Instance()->algo->fRadThick[i].GetRadThick(T.x, T.y), qp0, -1, wIn);
+ fit.L1AddMaterial(T, CbmL1::Instance()->algo->GetMaterialThickness(i, T.x, T.y), qp0, wIn);
+ fit.EnergyLossCorrection(T, CbmL1::Instance()->algo->GetMaterialThickness(i, T.x, T.y), qp0, -1, wIn);
}
else {
fit.L1AddMaterial(T, sta[i].materialInfo, qp0, wIn);
@@ -338,8 +338,8 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo)
fit.EnergyLossCorrection(T, fit.PipeRadThick, qp0, fvec(1.f), wIn);
}
if constexpr (L1Parameters::kIfUseRadLengthTable) {
- fit.L1AddMaterial(T, CbmL1::Instance()->algo->fRadThick[i].GetRadThick(T.x, T.y), qp0, wIn);
- fit.EnergyLossCorrection(T, CbmL1::Instance()->algo->fRadThick[i].GetRadThick(T.x, T.y), qp0, 1, wIn);
+ fit.L1AddMaterial(T, CbmL1::Instance()->algo->GetMaterialThickness(i, T.x, T.y), qp0, wIn);
+ fit.EnergyLossCorrection(T, CbmL1::Instance()->algo->GetMaterialThickness(i, T.x, T.y), qp0, 1, wIn);
}
else {
fit.L1AddMaterial(T, sta[i].materialInfo, qp0, wIn);
@@ -503,8 +503,8 @@ void CbmL1PFFitter::GetChiToVertex(vector<CbmStsTrack>& Tracks, vector<L1FieldRe
fit.L1AddPipeMaterial(T, T.qp, w);
fit.EnergyLossCorrection(T, fit.PipeRadThick, T.qp, fvec(1.f), w);
}
- fit.L1AddMaterial(T, CbmL1::Instance()->algo->fRadThick[iSt].GetRadThick(T.x, T.y), T.qp, w);
- fit.EnergyLossCorrection(T, CbmL1::Instance()->algo->fRadThick[iSt].GetRadThick(T.x, T.y), T.qp, fvec(1.f), w);
+ fit.L1AddMaterial(T, CbmL1::Instance()->algo->GetMaterialThickness(iSt, T.x, T.y), T.qp, w);
+ fit.EnergyLossCorrection(T, CbmL1::Instance()->algo->GetMaterialThickness(iSt, T.x, T.y), T.qp, fvec(1.f), w);
}
if (NMvdStations <= 0) {
fit.L1AddPipeMaterial(T, T.qp, ONE);