From 022b0112d7a57e33181a1ba5d800792189062f38 Mon Sep 17 00:00:00 2001
From: "se.gorbunov" <se.gorbunov@gsi.de>
Date: Wed, 9 Nov 2022 16:38:18 +0000
Subject: [PATCH] L1: remove option to use homogenious material
---
reco/L1/CMakeLists.txt | 2 +-
reco/L1/CbmL1.cxx | 38 ++---
reco/L1/CbmL1.h | 65 ---------
reco/L1/CbmL1Performance.cxx | 30 ++--
reco/L1/CbmL1ReadEvent.cxx | 16 +--
reco/L1/L1Algo/L1BaseStationInfo.cxx | 70 ++--------
reco/L1/L1Algo/L1BaseStationInfo.h | 68 +++------
reco/L1/L1Algo/L1BranchExtender.cxx | 41 +++---
reco/L1/L1Algo/L1CATrackFinder.cxx | 83 ++++-------
reco/L1/L1Algo/L1CloneMerger.cxx | 2 +-
reco/L1/L1Algo/L1Constants.h | 5 -
reco/L1/L1Algo/L1Fit.h | 7 +-
reco/L1/L1Algo/L1FitMaterial.cxx | 51 +------
reco/L1/L1Algo/L1InitManager.cxx | 11 +-
.../{L1MaterialInfo.cxx => L1Material.cxx} | 132 ++++++++++--------
.../L1Algo/{L1MaterialInfo.h => L1Material.h} | 42 +-----
reco/L1/L1Algo/L1Parameters.cxx | 4 +-
reco/L1/L1Algo/L1Parameters.h | 2 +-
reco/L1/L1Algo/L1Station.cxx | 8 +-
reco/L1/L1Algo/L1Station.h | 17 ++-
reco/L1/L1Algo/L1TrackFitter.cxx | 62 ++++----
reco/L1/L1Algo/L1TrackParFit.cxx | 27 ----
reco/L1/L1Algo/L1TrackParFit.h | 3 -
reco/L1/L1Algo/utils/L1AlgoDraw.cxx | 8 +-
reco/L1/ParticleFinder/CbmL1PFFitter.cxx | 28 ++--
25 files changed, 256 insertions(+), 566 deletions(-)
rename reco/L1/L1Algo/{L1MaterialInfo.cxx => L1Material.cxx} (66%)
rename reco/L1/L1Algo/{L1MaterialInfo.h => L1Material.h} (72%)
diff --git a/reco/L1/CMakeLists.txt b/reco/L1/CMakeLists.txt
index 7519440585..b4e445b438 100644
--- a/reco/L1/CMakeLists.txt
+++ b/reco/L1/CMakeLists.txt
@@ -52,7 +52,7 @@ set(SRCS
L1Algo/L1FitMaterial.cxx
L1Algo/L1Extrapolation.cxx
CbmL1MCTrack.cxx
- L1Algo/L1MaterialInfo.cxx
+ L1Algo/L1Material.cxx
L1Algo/L1UMeasurementInfo.cxx
L1Algo/L1XYMeasurementInfo.cxx
L1Algo/L1Field.cxx
diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx
index a734426621..ecabd503cc 100644
--- a/reco/L1/CbmL1.cxx
+++ b/reco/L1/CbmL1.cxx
@@ -464,9 +464,6 @@ InitStatus CbmL1::Init()
// *** MVD stations info ***
if (fUseMVD) {
auto materialTableMvd = ReadMaterialBudget(L1DetectorID::kMvd);
- auto correctionMvd = [this](L1Material& material, const L1MaterialInfo& homogenious) {
- this->ApplyCorrectionToMaterialMap<L1DetectorID::kMvd>(material, homogenious);
- };
for (int iSt = 0; iSt < fNMvdStationsGeom; ++iSt) {
auto stationInfo = L1BaseStationInfo(L1DetectorID::kMvd, iSt);
stationInfo.SetStationType(1); // MVD
@@ -478,8 +475,8 @@ InitStatus CbmL1::Init()
stationInfo.SetYmax(mvdInterface->GetYmax(iSt));
stationInfo.SetRmin(mvdInterface->GetRmin(iSt));
stationInfo.SetRmax(mvdInterface->GetRmax(iSt));
- stationInfo.SetMaterialSimple(mvdInterface->GetThickness(iSt), mvdInterface->GetRadLength(iSt));
- stationInfo.SetMaterialMap(std::move(materialTableMvd[iSt]), correctionMvd);
+ stationInfo.SetZthickness(mvdInterface->GetThickness(iSt));
+ stationInfo.SetMaterialMap(std::move(materialTableMvd[iSt]));
// TODO: The CA TF result is dependent from type of geometry settings. Should be understood (S.Zharko)
stationInfo.SetFrontBackStripsGeometry(
(fscal) mvdInterface->GetStripsStereoAngleFront(iSt), (fscal) mvdInterface->GetStripsSpatialRmsFront(iSt),
@@ -493,9 +490,6 @@ InitStatus CbmL1::Init()
// *** STS stations info ***
if (fUseSTS) {
auto materialTableSts = ReadMaterialBudget(L1DetectorID::kSts);
- auto correctionSts = [this](L1Material& material, const L1MaterialInfo& homogenious) {
- this->ApplyCorrectionToMaterialMap<L1DetectorID::kSts>(material, homogenious);
- };
for (int iSt = 0; iSt < fNStsStationsGeom; ++iSt) {
auto stationInfo = L1BaseStationInfo(L1DetectorID::kSts, iSt);
stationInfo.SetStationType(0); // STS
@@ -509,8 +503,8 @@ InitStatus CbmL1::Init()
stationInfo.SetYmax(stsInterface->GetYmax(iSt));
stationInfo.SetRmin(stsInterface->GetRmin(iSt));
stationInfo.SetRmax(stsInterface->GetRmax(iSt));
- stationInfo.SetMaterialSimple(stsInterface->GetThickness(iSt), stsInterface->GetRadLength(iSt));
- stationInfo.SetMaterialMap(std::move(materialTableSts[iSt]), correctionSts);
+ stationInfo.SetZthickness(stsInterface->GetThickness(iSt));
+ stationInfo.SetMaterialMap(std::move(materialTableSts[iSt]));
// TODO: The CA TF result is dependent from type of geometry settings. Should be understood (S.Zharko)
stationInfo.SetFrontBackStripsGeometry(
(fscal) stsInterface->GetStripsStereoAngleFront(iSt), (fscal) stsInterface->GetStripsSpatialRmsFront(iSt),
@@ -524,9 +518,6 @@ InitStatus CbmL1::Init()
// *** MuCh stations info ***
if (fUseMUCH) {
auto materialTableMuch = ReadMaterialBudget(L1DetectorID::kMuch);
- auto correctionMuch = [this](L1Material& material, const L1MaterialInfo& homogenious) {
- this->ApplyCorrectionToMaterialMap<L1DetectorID::kMuch>(material, homogenious);
- };
for (int iSt = 0; iSt < fNMuchStationsGeom; ++iSt) {
auto stationInfo = L1BaseStationInfo(L1DetectorID::kMuch, iSt);
stationInfo.SetStationType(2); // MuCh
@@ -540,8 +531,8 @@ InitStatus CbmL1::Init()
stationInfo.SetYmax(muchInterface->GetYmax(iSt));
stationInfo.SetRmin(muchInterface->GetRmin(iSt));
stationInfo.SetRmax(muchInterface->GetRmax(iSt));
- stationInfo.SetMaterialSimple(muchInterface->GetThickness(iSt), muchInterface->GetRadLength(iSt));
- stationInfo.SetMaterialMap(std::move(materialTableMuch[iSt]), correctionMuch);
+ stationInfo.SetZthickness(muchInterface->GetThickness(iSt));
+ stationInfo.SetMaterialMap(std::move(materialTableMuch[iSt]));
// TODO: The CA TF result is dependent from type of geometry settings. Should be understood (S.Zharko)
stationInfo.SetFrontBackStripsGeometry(
(fscal) muchInterface->GetStripsStereoAngleFront(iSt), (fscal) muchInterface->GetStripsSpatialRmsFront(iSt),
@@ -555,9 +546,6 @@ InitStatus CbmL1::Init()
// *** TRD stations info ***
if (fUseTRD) {
auto materialTableTrd = ReadMaterialBudget(L1DetectorID::kTrd);
- auto correctionTrd = [this](L1Material& material, const L1MaterialInfo& homogenious) {
- this->ApplyCorrectionToMaterialMap<L1DetectorID::kTrd>(material, homogenious);
- };
for (int iSt = 0; iSt < fNTrdStationsGeom; ++iSt) {
auto stationInfo = L1BaseStationInfo(L1DetectorID::kTrd, iSt);
stationInfo.SetStationType((iSt == 1 || iSt == 3) ? 6 : 3); // MuCh
@@ -571,8 +559,8 @@ InitStatus CbmL1::Init()
stationInfo.SetYmax(trdInterface->GetYmax(iSt));
stationInfo.SetRmin(trdInterface->GetRmin(iSt));
stationInfo.SetRmax(trdInterface->GetRmax(iSt));
- stationInfo.SetMaterialSimple(trdInterface->GetThickness(iSt), trdInterface->GetRadLength(iSt));
- stationInfo.SetMaterialMap(std::move(materialTableTrd[iSt]), correctionTrd);
+ stationInfo.SetZthickness(trdInterface->GetThickness(iSt));
+ stationInfo.SetMaterialMap(std::move(materialTableTrd[iSt]));
fscal trdFrontPhi = trdInterface->GetStripsStereoAngleFront(iSt);
fscal trdBackPhi = trdInterface->GetStripsStereoAngleBack(iSt);
fscal trdFrontSigma = trdInterface->GetStripsSpatialRmsFront(iSt);
@@ -596,9 +584,6 @@ InitStatus CbmL1::Init()
// *** TOF stations info ***
if (fUseTOF) {
auto materialTableTof = ReadMaterialBudget(L1DetectorID::kTof);
- auto correctionTof = [this](L1Material& material, const L1MaterialInfo& homogenious) {
- this->ApplyCorrectionToMaterialMap<L1DetectorID::kTof>(material, homogenious);
- };
for (int iSt = 0; iSt < fNTofStationsGeom; ++iSt) {
auto stationInfo = L1BaseStationInfo(L1DetectorID::kTof, iSt);
stationInfo.SetStationType(4);
@@ -609,9 +594,8 @@ InitStatus CbmL1::Init()
stationInfo.SetFieldStatus(0);
stationInfo.SetZ(tofInterface->GetZ(iSt));
auto thickness = tofInterface->GetThickness(iSt);
- auto radLength = tofInterface->GetRadLength(iSt);
- stationInfo.SetMaterialSimple(thickness, radLength);
- stationInfo.SetMaterialMap(std::move(materialTableTof[iSt]), correctionTof);
+ stationInfo.SetZthickness(thickness);
+ stationInfo.SetMaterialMap(std::move(materialTableTof[iSt]));
stationInfo.SetXmax(tofInterface->GetXmax(iSt));
stationInfo.SetYmax(tofInterface->GetYmax(iSt));
stationInfo.SetRmin(tofInterface->GetRmin(iSt));
@@ -1710,6 +1694,7 @@ std::vector<L1Material> CbmL1::ReadMaterialBudget(L1DetectorID detectorID)
result[iSt].SetRadThickBin(iBinX, iBinY, 0.01 * hStaRadLen->GetBinContent(iBinX, iBinY));
} // iBinX
} // iBinY
+ result[iSt].Repare();
LOG(info) << "- station " << iSt;
} // iSt
gFile = oldFile;
@@ -1721,6 +1706,7 @@ std::vector<L1Material> CbmL1::ReadMaterialBudget(L1DetectorID detectorID)
return result;
}
+
double CbmL1::boundedGaus(double sigma)
{
assert(sigma > 0. && std::isfinite(sigma));
diff --git a/reco/L1/CbmL1.h b/reco/L1/CbmL1.h
index bed97439cd..849f04022c 100644
--- a/reco/L1/CbmL1.h
+++ b/reco/L1/CbmL1.h
@@ -273,11 +273,6 @@ public:
/// Gets flag: to correct input hits on MC or not
bool GetCorrectHitsOnMC() const { return fIfCorrectHitsOnMC; }
- /// Correction function for the material budget map
- /// It fills bins with no statistics
- template<L1DetectorID detID>
- void ApplyCorrectionToMaterialMap(L1Material& material, const L1MaterialInfo& homogenious);
-
/// Utility to map the L1DetectorID items into detector names
static constexpr const char* GetDetectorName(L1DetectorID detectorID)
{
@@ -662,64 +657,4 @@ private:
};
-// ---------------------------------------------------------------------------------------------------------------------
-//
-template<L1DetectorID detID>
-void CbmL1::ApplyCorrectionToMaterialMap(L1Material& material, const L1MaterialInfo& homogenious)
-{
- // TODO: unify the correction function for all detectors
- [[maybe_unused]] float minVal = 0.;
- if constexpr (detID == L1DetectorID::kMuch) { minVal = 0.15f; }
- else if constexpr (detID == L1DetectorID::kTof || detID == L1DetectorID::kTrd) {
- minVal = 0.0015f;
- }
-
- // A bit ugly solution, but so can we avoid dependency on input maps file
- std::vector<float> keepRow {};
- if constexpr (detID != L1DetectorID::kSts) { keepRow.resize(material.GetNbins()); }
-
- for (int iBinX = 0; iBinX < material.GetNbins(); ++iBinX) {
- if constexpr (detID == L1DetectorID::kTof) { minVal = 0.0015f; }
- if constexpr (detID != L1DetectorID::kSts) {
- for (int iBinY = 0; iBinY < material.GetNbins(); ++iBinY) {
- keepRow[iBinY] = material.GetRadThickBin(iBinX, iBinY);
- }
- }
- for (int iBinY = 0; iBinY < material.GetNbins(); ++iBinY) {
- if constexpr (detID == L1DetectorID::kMvd) {
- // Correction for holes in the material map
- if (material.GetRadThickBin(iBinX, iBinY) < homogenious.RadThick[0]) {
- if (iBinY > 0 && iBinY < material.GetNbins() - 1) {
- material.SetRadThickBin(iBinX, iBinY, TMath::Min(keepRow[iBinY - 1], keepRow[iBinY + 1]));
- }
- }
-
- // Correction for the hard-coded value of RadThick of MVD stations
- if (material.GetRadThickBin(iBinX, iBinY) < 0.0015) { material.SetRadThickBin(iBinX, iBinY, 0.0015); }
- }
- else if constexpr (detID == L1DetectorID::kSts) {
- if (material.GetRadThickBin(iBinX, iBinY) < homogenious.RadThick[0]) {
- material.SetRadThickBin(iBinX, iBinY, homogenious.RadThick[0]);
- }
- }
- else if constexpr (detID == L1DetectorID::kMuch || detID == L1DetectorID::kTrd || detID == L1DetectorID::kTof) {
- // Correction for holes in the material map
- if (L1Algo::TrackingMode::kGlobal != fTrackingMode) {
- if ((iBinY > 0) && (iBinY < material.GetNbins() - 1)) {
- material.SetRadThickBin(iBinX, iBinY, TMath::Min(keepRow[iBinY - 1], keepRow[iBinY + 1]));
- }
- }
- float val = material.GetRadThickBin(iBinX, iBinY);
- if (val > 0.0015) { // remember last non-zero value
- minVal = val;
- }
- else { // empty bin with no statistics, fill it with the neighbours value
- material.SetRadThickBin(iBinX, iBinY, minVal);
- }
- }
- }
- }
-}
-
-
#endif //_CbmL1_h_
diff --git a/reco/L1/CbmL1Performance.cxx b/reco/L1/CbmL1Performance.cxx
index 5b855ce4a4..da1c6784a8 100644
--- a/reco/L1/CbmL1Performance.cxx
+++ b/reco/L1/CbmL1Performance.cxx
@@ -823,8 +823,8 @@ void CbmL1::HistoPerformance() // TODO: check if works correctly. Change vHitFa
CbmL1HitStore& h2 = fvHitStore[prtra->Hits[1]];
h_ghost_fstation->Fill(h1.iStation);
h_ghost_r->Fill(sqrt(fabs(h1.x * h1.x + h1.y * h1.y)));
- double z1 = fpAlgo->GetParameters()->GetStation(h1.iStation).z[0];
- double z2 = fpAlgo->GetParameters()->GetStation(h2.iStation).z[0];
+ double z1 = fpAlgo->GetParameters()->GetStation(h1.iStation).fZ[0];
+ double z2 = fpAlgo->GetParameters()->GetStation(h2.iStation).fZ[0];
if (fabs(z2 - z1) > 1.e-4) {
h_ghost_tx->Fill((h2.x - h1.x) / (z2 - z1));
h_ghost_ty->Fill((h2.y - h1.y) / (z2 - z1));
@@ -1025,8 +1025,8 @@ void CbmL1::HistoPerformance() // TODO: check if works correctly. Change vHitFa
// CbmL1HitStore &ph21 = fvHitStore[mtra.Hits[0]];
// CbmL1HitStore &ph22 = fvHitStore[mtra.Hits[1]];
- // double z21 = fpAlgo->GetParameters()->GetStation(ph21.iStation).z[0];
- // double z22 = fpAlgo->GetParameters()->GetStation(ph22.iStation).z[0];
+ // double z21 = fpAlgo->GetParameters()->GetStation(ph21.iStation).fZ[0];
+ // double z22 = fpAlgo->GetParameters()->GetStation(ph22.iStation).fZ[0];
// if( fabs(z22-z21)>1.e-4 ){
// h_notfound_tx->Fill((ph22.x-ph21.x)/(z22-z21));
// h_notfound_ty->Fill((ph22.y-ph21.y)/(z22-z21));
@@ -1395,14 +1395,16 @@ void CbmL1::TrackFitPerformance()
L1Extrapolate(trPar, mc.z, trPar.qp, fld);
// add material
const int fSta = fvHitStore[it->Hits[0]].iStation;
- const int dir = int((mc.z - fpAlgo->GetParameters()->GetStation(fSta).z[0])
- / fabs(mc.z - fpAlgo->GetParameters()->GetStation(fSta).z[0]));
- // if (abs(mc.z - fpAlgo->GetParameters()->GetStation(fSta).z[0]) > 10.) continue; // can't extrapolate on large distance
+ const int dir = int((mc.z - fpAlgo->GetParameters()->GetStation(fSta).fZ[0])
+ / fabs(mc.z - fpAlgo->GetParameters()->GetStation(fSta).fZ[0]));
+ // if (abs(mc.z - fpAlgo->GetParameters()->GetStation(fSta).fZ[0]) > 10.) continue; // can't extrapolate on large distance
for (int iSta = fSta /*+dir*/; (iSta >= 0) && (iSta < fNStations)
- && (dir * (mc.z - fpAlgo->GetParameters()->GetStation(iSta).z[0]) > 0);
+ && (dir * (mc.z - fpAlgo->GetParameters()->GetStation(iSta).fZ[0]) > 0);
iSta += dir) {
// cout << iSta << " " << dir << endl;
- fit.L1AddMaterial(trPar, fpAlgo->GetParameters()->GetStation(iSta).materialInfo, trPar.qp, 1);
+ fit.L1AddMaterial(trPar, fpAlgo->GetParameters()->GetMaterialThickness(iSta, trPar.x, trPar.y), trPar.qp,
+ fvec::One());
+
if (iSta + dir == fNMvdStations - 1) fit.L1AddPipeMaterial(trPar, trPar.qp, 1);
}
}
@@ -1448,15 +1450,15 @@ void CbmL1::TrackFitPerformance()
// add material
const int fSta = fvHitStore[it->Hits[0]].iStation;
- const int dir = (mc.z - fpAlgo->GetParameters()->GetStation(fSta).z[0])
- / abs(mc.z - fpAlgo->GetParameters()->GetStation(fSta).z[0]);
- // if (abs(mc.z - fpAlgo->GetParameters()->GetStation(fSta].z[0]) > 10.) continue; // can't extrapolate on large distance
+ const int dir = (mc.z - fpAlgo->GetParameters()->GetStation(fSta).fZ[0])
+ / abs(mc.z - fpAlgo->GetParameters()->GetStation(fSta).fZ[0]);
+ // if (abs(mc.z - fpAlgo->GetParameters()->GetStation(fSta].fZ[0]) > 10.) continue; // can't extrapolate on large distance
for (int iSta = fSta + dir; (iSta >= 0) && (iSta < fNStations)
- && (dir * (mc.z - fpAlgo->GetParameters()->GetStation(iSta).z[0]) > 0);
+ && (dir * (mc.z - fpAlgo->GetParameters()->GetStation(iSta).fZ[0]) > 0);
iSta += dir) {
- L1Extrapolate(trPar, fpAlgo->GetParameters()->GetStation(iSta).z[0], trPar.qp, fld);
+ L1Extrapolate(trPar, fpAlgo->GetParameters()->GetStation(iSta).fZ[0], trPar.qp, fld);
fit.L1AddMaterial(trPar, fpAlgo->GetParameters()->GetMaterialThickness(iSta, trPar.x, trPar.y), trPar.qp,
1);
fit.EnergyLossCorrection(trPar, fpAlgo->GetParameters()->GetMaterialThickness(iSta, trPar.x, trPar.y),
diff --git a/reco/L1/CbmL1ReadEvent.cxx b/reco/L1/CbmL1ReadEvent.cxx
index c4695c66df..79493eeafc 100644
--- a/reco/L1/CbmL1ReadEvent.cxx
+++ b/reco/L1/CbmL1ReadEvent.cxx
@@ -439,7 +439,7 @@ void CbmL1::ReadEvent(float& TsStart, float& TsLength, float& /*TsOverlap*/, int
// due to a problem in transport
int iStActive = fpAlgo->GetParameters()->GetStationIndexActive(iSt, L1DetectorID::kMvd);
if (iStActive == -1) { continue; }
- double d = (MC.zIn - sta[iStActive].z[0]);
+ double d = (MC.zIn - sta[iStActive].fZ[0]);
if (fabs(d) < fabs(bestDist)) {
bestDist = d;
MC.iStation = iStActive;
@@ -480,7 +480,7 @@ void CbmL1::ReadEvent(float& TsStart, float& TsLength, float& /*TsOverlap*/, int
if (iStActive == -1) { continue; }
// use z_in since z_out is sometimes very wrong
// due to a problem in transport
- double d = (MC.zIn - sta[iStActive].z[0]);
+ double d = (MC.zIn - sta[iStActive].fZ[0]);
if (fabs(d) < fabs(bestDist)) {
bestDist = d;
MC.iStation = iStActive;
@@ -514,7 +514,7 @@ void CbmL1::ReadEvent(float& TsStart, float& TsLength, float& /*TsOverlap*/, int
for (Int_t iSt = 0; iSt < fNMuchStationsGeom; iSt++) {
int iStActive = fpAlgo->GetParameters()->GetStationIndexActive(iSt, L1DetectorID::kMuch);
assert(iStActive != -1);
- if (MC.z > sta[iStActive].z[0] - 2.5) { MC.iStation = iStActive; }
+ if (MC.z > sta[iStActive].fZ[0] - 2.5) { MC.iStation = iStActive; }
}
assert(MC.iStation >= 0);
auto itTrack = fmMCTracksLinksMap.find(CbmL1LinkKey(MC.ID, iEvent, iFile));
@@ -541,7 +541,7 @@ void CbmL1::ReadEvent(float& TsStart, float& TsLength, float& /*TsOverlap*/, int
int iStActive = fpAlgo->GetParameters()->GetStationIndexActive(iSt, L1DetectorID::kTrd);
assert(iStActive != -1);
- if (MC.z > sta[iStActive].z[0] - 20.0) { MC.iStation = iStActive; }
+ if (MC.z > sta[iStActive].fZ[0] - 20.0) { MC.iStation = iStActive; }
}
assert(MC.iStation >= 0);
auto itTrack = fmMCTracksLinksMap.find(CbmL1LinkKey(MC.ID, iEvent, iFile));
@@ -608,17 +608,17 @@ void CbmL1::ReadEvent(float& TsStart, float& TsLength, float& /*TsOverlap*/, int
for (int iSt = 0; iSt < fNTofStationsGeom; iSt++) {
int iStActive = fpAlgo->GetParameters()->GetStationIndexActive(iSt, L1DetectorID::kTof);
if (iStActive == -1) { continue; }
- if (fabs(MC.z - sta[iStActive].z[0]) < dist) {
- dist = fabs(MC.z - sta[iStActive].z[0]);
+ if (fabs(MC.z - sta[iStActive].fZ[0]) < dist) {
+ dist = fabs(MC.z - sta[iStActive].fZ[0]);
iSta = iSt;
}
}
MC.iStation = fpAlgo->GetParameters()->GetStationIndexActive(iSta, L1DetectorID::kTof);
assert(MC.iStation >= 0);
if (iSta >= 0) {
- if (fabs(sta[iSta].z[0] - MC.z) < TofPointToTrackdZ[iSta][iTrack]) {
+ if (fabs(sta[iSta].fZ[0] - MC.z) < TofPointToTrackdZ[iSta][iTrack]) {
fTofPointToTrack[iSta][iTrack] = iMC;
- TofPointToTrackdZ[iSta][iTrack] = fabs(sta[iSta].z[0] - MC.z);
+ TofPointToTrackdZ[iSta][iTrack] = fabs(sta[iSta].fZ[0] - MC.z);
}
}
}
diff --git a/reco/L1/L1Algo/L1BaseStationInfo.cxx b/reco/L1/L1Algo/L1BaseStationInfo.cxx
index 6165bc7817..ba2c03518e 100644
--- a/reco/L1/L1Algo/L1BaseStationInfo.cxx
+++ b/reco/L1/L1Algo/L1BaseStationInfo.cxx
@@ -335,25 +335,21 @@ void L1BaseStationInfo::SetFrontBackStripsGeometry(double frontPhi, double front
//------------------------------------------------------------------------------------------------------------------------
//
-void L1BaseStationInfo::SetMaterialSimple(double inThickness, double inRL)
+void L1BaseStationInfo::SetZthickness(double inThickness)
{
//L1MASSERT(0, inRL, "Attempt of entering zero inRL (radiational length) value");
- fL1Station.materialInfo.thick = inThickness;
- fL1Station.materialInfo.RL = inRL;
- fL1Station.materialInfo.RadThick = fL1Station.materialInfo.thick / fL1Station.materialInfo.RL;
- fL1Station.materialInfo.logRadThick = log(fL1Station.materialInfo.RadThick);
- fInitController.SetFlag(EInitKey::kMaterialInfo);
+ fL1Station.fZthick = inThickness;
+ fInitController.SetFlag(EInitKey::kZthickness);
}
//------------------------------------------------------------------------------------------------------------------------
//
-void L1BaseStationInfo::SetMaterialMap(const L1Material& thicknessMap,
- const std::function<void(L1Material& thicknessMap)>& correction)
+void L1BaseStationInfo::SetMaterialMap(const L1Material& thicknessMap)
{
if (!fInitController.GetFlag(EInitKey::kThicknessMap)) {
fThicknessMap = thicknessMap;
- if (correction) { correction(fThicknessMap); }
+ fThicknessMap.Repare();
fInitController.SetFlag(EInitKey::kThicknessMap);
}
else {
@@ -361,39 +357,13 @@ void L1BaseStationInfo::SetMaterialMap(const L1Material& thicknessMap,
}
}
-//------------------------------------------------------------------------------------------------------------------------
-//
-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
+void L1BaseStationInfo::SetMaterialMap(L1Material&& thicknessMap) noexcept
{
if (!fInitController.GetFlag(EInitKey::kThicknessMap)) {
fThicknessMap = std::move(thicknessMap);
- if (correction) { correction(fThicknessMap); }
fInitController.SetFlag(EInitKey::kThicknessMap);
}
else {
@@ -401,30 +371,6 @@ void L1BaseStationInfo::SetMaterialMap(L1Material&& thicknessMap,
}
}
-//------------------------------------------------------------------------------------------------------------------------
-//
-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";
- }
-}
//------------------------------------------------------------------------------------------------------------------------
//
@@ -496,8 +442,8 @@ void L1BaseStationInfo::SetTrackingStatus(bool flag)
//
void L1BaseStationInfo::SetZ(double inZ)
{
- fL1Station.z = inZ; // setting simd vector of single-precision floats, which is passed to high performanced L1Algo
- fZPos = inZ; // setting precised value to use in field approximation etc
+ fL1Station.fZ = inZ; // setting simd vector of single-precision floats, which is passed to high performanced L1Algo
+ fZPos = inZ; // setting precised value to use in field approximation etc
fInitController.SetFlag(EInitKey::kZ);
}
diff --git a/reco/L1/L1Algo/L1BaseStationInfo.h b/reco/L1/L1Algo/L1BaseStationInfo.h
index 45910b0ec3..0f3a51fad3 100644
--- a/reco/L1/L1Algo/L1BaseStationInfo.h
+++ b/reco/L1/L1Algo/L1BaseStationInfo.h
@@ -46,20 +46,20 @@ public:
kXmax, ///< max size in X direction
kYmax, ///< max size in Y direction
// L1Station initialization
- kType, ///< station type
- kTimeInfo, ///< if time info is used (flag)
- kFieldStatus, ///< if station is placed in field (flag)
- kZ, ///< z coordinate of the station position
- kRmin, ///< internal radius of station (gap size)
- kRmax, ///< exteranl radius of 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, ///<
- kStripsBackPhi, ///<
- kStripsBackSigma, ///<
- kTimeResolution, ///< time resolution
+ kType, ///< station type
+ kTimeInfo, ///< if time info is used (flag)
+ kFieldStatus, ///< if station is placed in field (flag)
+ kZ, ///< z coordinate of the station position
+ kRmin, ///< internal radius of station (gap size)
+ kRmax, ///< exteranl radius of station
+ kZthickness, ///< Z thickness of the station
+ kThicknessMap, ///< thickness map of the station (optional?)
+ kFieldSlice, ///< L1Station.L1FieldSlice object initialization
+ kStripsFrontPhi, ///< strips geometry initialization
+ kStripsFrontSigma, ///<
+ kStripsBackPhi, ///<
+ kStripsBackSigma, ///<
+ kTimeResolution, ///< time resolution
// The last item is equal to the number of bits in fInitFlags
kEnd
};
@@ -131,16 +131,7 @@ public:
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 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; }
+ fvec GetZthickness() const { return fL1Station.fZthick; }
/// Gets min transverse size of the station [cm]
fvec GetRmin() const { return fL1Station.Rmin; }
@@ -170,7 +161,7 @@ public:
double GetZdouble() const { return fZPos; }
/// Gets SIMD vectorized z position of the station [cm]
- fvec GetZsimdVec() const { return fL1Station.z; }
+ fvec GetZsimdVec() const { return fL1Station.fZ; }
/// Prints registered fields
/// verbosity = 0: print only station id, detector id and address in one line
@@ -201,36 +192,15 @@ public:
/// 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);
+ void SetZthickness(double thickness);
/// 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;
+ void SetMaterialMap(const L1Material& thicknessMap);
/// 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;
+ void SetMaterialMap(L1Material&& thicknessMap) noexcept;
/// Sets max transverse size of the station []
void SetRmax(double inRmax);
diff --git a/reco/L1/L1Algo/L1BranchExtender.cxx b/reco/L1/L1Algo/L1BranchExtender.cxx
index 3ed23c7194..2c3d46bf79 100644
--- a/reco/L1/L1Algo/L1BranchExtender.cxx
+++ b/reco/L1/L1Algo/L1BranchExtender.cxx
@@ -100,9 +100,9 @@ void L1Algo::BranchFitterFast(const L1Branch& t, L1TrackPar& T, const bool dir,
L1FieldValue fldB0, fldB1, fldB2 _fvecalignment;
L1FieldRegion fld _fvecalignment;
- fvec fldZ0 = sta1.z; // suppose field is smoth
- fvec fldZ1 = sta2.z;
- fvec fldZ2 = sta0.z;
+ fvec fldZ0 = sta1.fZ; // suppose field is smoth
+ fvec fldZ1 = sta2.fZ;
+ fvec fldZ2 = sta0.fZ;
sta1.fieldSlice.GetFieldValue(x1, y1, fldB0);
@@ -111,13 +111,13 @@ void L1Algo::BranchFitterFast(const L1Branch& t, L1TrackPar& T, const bool dir,
fld.Set(fldB2, fldZ2, fldB1, fldZ1, fldB0, fldZ0);
- int ista_prev = ista1;
- int ista = ista2;
+ //int ista_prev = ista1;
+ int ista = ista2;
for (int i = iFirstHit + step; step * i <= step * iLastHit; i += step) {
const L1Hit& hit = fInputData.GetHit(hits[i]);
- ista_prev = ista;
- ista = hit.iSt;
+ //ista_prev = ista;
+ ista = hit.iSt;
const L1Station& sta = fParameters.GetStation(ista);
@@ -130,11 +130,10 @@ void L1Algo::BranchFitterFast(const L1Branch& t, L1TrackPar& T, const bool dir,
L1Extrapolate(T, hit.z, qp0, fld);
}
L1ExtrapolateTime(T, dz);
-
- fit.L1AddMaterial(T, sta.materialInfo, qp0, 1);
- if ((step * ista <= step * (fNstationsBeforePipe + (step + 1) / 2 - 1))
- && (step * ista_prev >= step * (fNstationsBeforePipe + (step + 1) / 2 - 1 - step)))
- fit.L1AddPipeMaterial(T, qp0, 1);
+ fit.L1AddMaterial(T, fParameters.GetMaterialThickness(ista, T.x, T.y), qp0, fvec::One());
+ //if ((step * ista <= step * (fNstationsBeforePipe + (step + 1) / 2 - 1))
+ // && (step * ista_prev >= step * (fNstationsBeforePipe + (step + 1) / 2 - 1 - step)))
+ //fit.L1AddPipeMaterial(T, qp0, 1);
fvec u = hit.u;
fvec v = hit.v;
@@ -158,7 +157,7 @@ void L1Algo::BranchFitterFast(const L1Branch& t, L1TrackPar& T, const bool dir,
auto [x, y] = sta.ConvUVtoXY<fvec>(u, v);
sta.fieldSlice.GetFieldValue(x, y, fldB2);
- fldZ2 = sta.z;
+ fldZ2 = sta.fZ;
fld.Set(fldB2, fldZ2, fldB1, fldZ1, fldB0, fldZ0);
} // i
@@ -218,9 +217,9 @@ void L1Algo::FindMoreHits(L1Branch& t, L1TrackPar& T, const bool dir,
L1FieldValue fldB0, fldB1, fldB2 _fvecalignment;
L1FieldRegion fld _fvecalignment;
- fvec fldZ0 = sta1.z;
- fvec fldZ1 = sta2.z;
- fvec fldZ2 = sta0.z;
+ fvec fldZ0 = sta1.fZ;
+ fvec fldZ1 = sta2.fZ;
+ fvec fldZ2 = sta0.fZ;
sta1.fieldSlice.GetFieldValue(x1, y1, fldB0);
sta2.fieldSlice.GetFieldValue(x2, y2, fldB1);
@@ -238,11 +237,11 @@ void L1Algo::FindMoreHits(L1Branch& t, L1TrackPar& T, const bool dir,
const L1Station& sta = fParameters.GetStation(ista);
- fvec dz = sta.z - T.z;
+ fvec dz = sta.fZ - T.z;
- if (kMcbm == fTrackingMode) { L1ExtrapolateLine(T, sta.z); }
+ if (kMcbm == fTrackingMode) { L1ExtrapolateLine(T, sta.fZ); }
else {
- L1Extrapolate(T, sta.z, qp0, fld);
+ L1Extrapolate(T, sta.fZ, qp0, fld);
}
L1ExtrapolateTime(T, dz);
@@ -316,7 +315,7 @@ void L1Algo::FindMoreHits(L1Branch& t, L1TrackPar& T, const bool dir,
L1ExtrapolateTime(T, dz1);
L1ExtrapolateLine(T, z);
- fit.L1AddMaterial(T, sta.materialInfo, qp0, 1);
+ fit.L1AddMaterial(T, fParameters.GetMaterialThickness(ista, T.x, T.y), qp0, fvec::One());
L1UMeasurementInfo info = sta.frontInfo;
@@ -335,7 +334,7 @@ void L1Algo::FindMoreHits(L1Branch& t, L1TrackPar& T, const bool dir,
fldZ0 = fldZ1;
fldZ1 = fldZ2;
sta.fieldSlice.GetFieldValue(x, y, fldB2);
- fldZ2 = sta.z;
+ fldZ2 = sta.fZ;
fld.Set(fldB2, fldZ2, fldB1, fldZ1, fldB0, fldZ0);
}
diff --git a/reco/L1/L1Algo/L1CATrackFinder.cxx b/reco/L1/L1Algo/L1CATrackFinder.cxx
index 4cd4ac7a5a..2a8c9e8b68 100644
--- a/reco/L1/L1Algo/L1CATrackFinder.cxx
+++ b/reco/L1/L1Algo/L1CATrackFinder.cxx
@@ -251,14 +251,14 @@ inline void L1Algo::findSingletsStep1( /// input 1st stage of singlet search
const fvec ty = (yl - fTargY) * dzli;
L1FieldValue b00, b01, b10, b11, b12;
- fld0Sta0.fieldSlice.GetFieldValue(fTargX + tx * (fld0Sta0.z - fTargZ), fTargY + ty * (fld0Sta0.z - fTargZ), b00);
- fld0Sta1.fieldSlice.GetFieldValue(fTargX + tx * (fld0Sta1.z - fTargZ), fTargY + ty * (fld0Sta1.z - fTargZ), b01);
- fld1Sta0.fieldSlice.GetFieldValue(fTargX + tx * (fld1Sta0.z - fTargZ), fTargY + ty * (fld1Sta0.z - fTargZ), b10);
- fld1Sta1.fieldSlice.GetFieldValue(fTargX + tx * (fld1Sta1.z - fTargZ), fTargY + ty * (fld1Sta1.z - fTargZ), b11);
- fld1Sta2.fieldSlice.GetFieldValue(fTargX + tx * (fld1Sta2.z - fTargZ), fTargY + ty * (fld1Sta2.z - fTargZ), b12);
+ fld0Sta0.fieldSlice.GetFieldValue(fTargX + tx * (fld0Sta0.fZ - fTargZ), fTargY + ty * (fld0Sta0.fZ - fTargZ), b00);
+ fld0Sta1.fieldSlice.GetFieldValue(fTargX + tx * (fld0Sta1.fZ - fTargZ), fTargY + ty * (fld0Sta1.fZ - fTargZ), b01);
+ fld1Sta0.fieldSlice.GetFieldValue(fTargX + tx * (fld1Sta0.fZ - fTargZ), fTargY + ty * (fld1Sta0.fZ - fTargZ), b10);
+ fld1Sta1.fieldSlice.GetFieldValue(fTargX + tx * (fld1Sta1.fZ - fTargZ), fTargY + ty * (fld1Sta1.fZ - fTargZ), b11);
+ fld1Sta2.fieldSlice.GetFieldValue(fTargX + tx * (fld1Sta2.fZ - fTargZ), fTargY + ty * (fld1Sta2.fZ - fTargZ), b12);
- fld0.Set(fTargB, fTargZ, b00, fld0Sta0.z, b01, fld0Sta1.z);
- fld1.Set(b10, fld1Sta0.z, b11, fld1Sta1.z, b12, fld1Sta2.z);
+ fld0.Set(fTargB, fTargZ, b00, fld0Sta0.fZ, b01, fld0Sta1.fZ);
+ fld1.Set(b10, fld1Sta0.fZ, b11, fld1Sta1.fZ, b12, fld1Sta2.fZ);
T.chi2 = 0.;
T.NDF = (fpCurrentIteration->GetPrimaryFlag()) ? fvec(2.) : fvec(0.);
@@ -338,34 +338,27 @@ inline void L1Algo::findSingletsStep1( /// input 1st stage of singlet search
//assert(T.IsConsistent(true, -1));
- if constexpr (L1Constants::control::kIfUseRadLengthTable) {
- if (kMcbm == fTrackingMode) {
- fit.L1AddThickMaterial(T, fParameters.GetMaterialThickness(istal, T.x, T.y), fMaxInvMom, fvec::One(),
- stal.materialInfo.thick, 1);
- }
- else if (kGlobal == fTrackingMode) {
- fit.L1AddMaterial(T, fParameters.GetMaterialThickness(istal, T.x, T.y), fMaxInvMom, fvec::One());
- }
- else {
- fit.L1AddMaterial(T, fParameters.GetMaterialThickness(istal, T.x, T.y), fMaxInvMom, fvec::One());
- }
+ if (kMcbm == fTrackingMode) {
+ fit.L1AddThickMaterial(T, fParameters.GetMaterialThickness(istal, T.x, T.y), fMaxInvMom, fvec::One(),
+ stal.fZthick, 1);
}
else {
- fit.L1AddMaterial(T, stal.materialInfo, fMaxInvMom, fvec::One());
+ fit.L1AddMaterial(T, fParameters.GetMaterialThickness(istal, T.x, T.y), fMaxInvMom, fvec::One());
}
+
//if ((istam >= fNstationsBeforePipe) && (istal <= fNstationsBeforePipe - 1)) {
//fit.L1AddPipeMaterial(T, fMaxInvMom, fvec::One());
//}
//assert(T.IsConsistent(true, -1));
- fvec dz = stam.z - zl;
+ fvec dz = stam.fZ - zl;
L1ExtrapolateTime(T, dz, stam.timeInfo);
// extrapolate to the middle hit
- if (istam < fNfieldStations) { L1Extrapolate0(T, stam.z, fld0); }
+ if (istam < fNfieldStations) { L1Extrapolate0(T, stam.fZ, fld0); }
else {
- L1ExtrapolateLine(T, stam.z); // TODO: fld1 doesn't work!
+ L1ExtrapolateLine(T, stam.fZ); // TODO: fld1 doesn't work!
}
// assert(T.IsConsistent(true, -1));
@@ -689,26 +682,22 @@ inline void L1Algo::findTripletsStep0( // input
// assert(T2.IsConsistent(true, n2_4));
- if constexpr (L1Constants::control::kIfUseRadLengthTable) {
- if (kMcbm == fTrackingMode) {
- fit.L1AddThickMaterial(T2, fParameters.GetMaterialThickness(istam, T2.x, T2.y), fMaxInvMom, fvec::One(),
- stam.materialInfo.thick, 1);
- }
- else if (kGlobal == fTrackingMode) {
- fit.L1AddMaterial(T2, fParameters.GetMaterialThickness(istam, T2.x, T2.y), fMaxInvMom, fvec::One());
- }
- else {
- fit.L1AddMaterial(T2, fParameters.GetMaterialThickness(istam, T2.x, T2.y), T2.qp, fvec::One());
- }
+ if (kMcbm == fTrackingMode) {
+ fit.L1AddThickMaterial(T2, fParameters.GetMaterialThickness(istam, T2.x, T2.y), fMaxInvMom, fvec::One(),
+ stam.fZthick, 1);
+ }
+ else if (kGlobal == fTrackingMode) {
+ fit.L1AddMaterial(T2, fParameters.GetMaterialThickness(istam, T2.x, T2.y), fMaxInvMom, fvec::One());
}
else {
- fit.L1AddMaterial(T2, stam.materialInfo, T2.qp, fvec::One());
+ fit.L1AddMaterial(T2, fParameters.GetMaterialThickness(istam, T2.x, T2.y), T2.qp, fvec::One());
}
+
//if ((istar >= fNstationsBeforePipe) && (istam <= fNstationsBeforePipe - 1)) { fit.L1AddPipeMaterial(T2, T2.qp, 1); }
// assert(T2.IsConsistent(true, n2_4));
- fvec dz2 = star.z - T2.z;
+ fvec dz2 = star.fZ - T2.z;
L1ExtrapolateTime(T2, dz2, stam.timeInfo);
// assert(T2.IsConsistent(true, n2_4));
@@ -716,10 +705,10 @@ inline void L1Algo::findTripletsStep0( // input
// extrapolate to the right hit station
if (istar <= fNfieldStations) {
- L1Extrapolate(T2, star.z, T2.qp, f2); // Full extrapolation in the magnetic field
+ L1Extrapolate(T2, star.fZ, T2.qp, f2); // Full extrapolation in the magnetic field
}
else {
- L1ExtrapolateLine(T2, star.z); // Extrapolation with line ()
+ L1ExtrapolateLine(T2, star.fZ); // Extrapolation with line ()
}
// assert(T2.IsConsistent(true, n2_4));
@@ -1039,12 +1028,12 @@ inline void L1Algo::findTripletsStep2(Tindex n3, int istal, int istam, int istar
fvec tx[3] = {(x[1] - x[0]) / (z[1] - z[0]), (x[2] - x[0]) / (z[2] - z[0]), (x[2] - x[1]) / (z[2] - z[1])};
fvec ty[3] = {(y[1] - y[0]) / (z[1] - z[0]), (y[2] - y[0]) / (z[2] - z[0]), (y[2] - y[1]) / (z[2] - z[1])};
for (int ih = 0; ih < NHits; ++ih) {
- fvec dz = (sta[ih].z - z[ih]);
+ fvec dz = (sta[ih].fZ - z[ih]);
sta[ih].fieldSlice.GetFieldValue(x[ih] + tx[ih] * dz, y[ih] + ty[ih] * dz, B[ih]);
};
- fld.Set(B[0], sta[0].z, B[1], sta[1].z, B[2], sta[2].z);
- fldTarget.Set(fTargB, fTargZ, B[0], sta[0].z, B[1], sta[1].z);
+ fld.Set(B[0], sta[0].fZ, B[1], sta[1].fZ, B[2], sta[2].fZ);
+ fldTarget.Set(fTargB, fTargZ, B[0], sta[0].fZ, B[1], sta[1].fZ);
L1TrackPar& T = fit.fTr;
@@ -1097,12 +1086,7 @@ inline void L1Algo::findTripletsStep2(Tindex n3, int istal, int istam, int istar
for (int ih = 1; ih < NHits; ++ih) {
fit.Extrapolate(z[ih], qp0, fld, fvec::One());
- if constexpr (L1Constants::control::kIfUseRadLengthTable) {
- fit.AddMaterial(fParameters.GetMaterialThickness(ista[ih], T.x, T.y), qp0, fvec::One());
- }
- else {
- fit.AddMaterial(sta[ih].materialInfo, qp0, fvec::One());
- }
+ fit.AddMaterial(fParameters.GetMaterialThickness(ista[ih], T.x, T.y), qp0, fvec::One());
//if (ista[ih] == fNstationsBeforePipe) { fit.AddPipeMaterial(qp0, fvec::One()); }
fit.Filter(sta[ih].frontInfo, u[ih], du2[ih], fvec::One());
fit.Filter(sta[ih].backInfo, v[ih], dv2[ih], fvec::One());
@@ -1135,12 +1119,7 @@ inline void L1Algo::findTripletsStep2(Tindex n3, int istal, int istam, int istar
for (int ih = NHits - 2; ih >= 0; --ih) {
fit.Extrapolate(z[ih], qp0, fld, fvec::One());
- if constexpr (L1Constants::control::kIfUseRadLengthTable) {
- fit.AddMaterial(fParameters.GetMaterialThickness(ista[ih], T.x, T.y), qp0, fvec::One());
- }
- else {
- fit.AddMaterial(sta[ih].materialInfo, qp0, fvec::One());
- }
+ fit.AddMaterial(fParameters.GetMaterialThickness(ista[ih], T.x, T.y), qp0, fvec::One());
//if (ista[ih] == fNstationsBeforePipe - 1) { fit.AddPipeMaterial(qp0, fvec::One()); }
fit.Filter(sta[ih].frontInfo, u[ih], du2[ih], fvec::One());
fit.Filter(sta[ih].backInfo, v[ih], dv2[ih], fvec::One());
diff --git a/reco/L1/L1Algo/L1CloneMerger.cxx b/reco/L1/L1Algo/L1CloneMerger.cxx
index 33d239964a..9d2d28b4ad 100644
--- a/reco/L1/L1Algo/L1CloneMerger.cxx
+++ b/reco/L1/L1Algo/L1CloneMerger.cxx
@@ -169,7 +169,7 @@ void L1CloneMerger::Exec(L1Vector<L1Track>& extTracks, L1Vector<L1HitIndex_t>& e
else
stam = staf - 1;
- fvec zm = frAlgo.GetParameters()->GetStation(stam).z;
+ fvec zm = frAlgo.GetParameters()->GetStation(stam).fZ;
fvec xm = fvec(0.5) * (Tf.x + Tf.tx * (zm - Tf.z) + Tb.x + Tb.tx * (zm - Tb.z));
fvec ym = fvec(0.5) * (Tf.y + Tf.ty * (zm - Tf.z) + Tb.y + Tb.ty * (zm - Tb.z));
frAlgo.GetParameters()->GetStation(stam).fieldSlice.GetFieldValue(xm, ym, fBm);
diff --git a/reco/L1/L1Algo/L1Constants.h b/reco/L1/L1Algo/L1Constants.h
index fa687d72e5..77d3a4c03c 100644
--- a/reco/L1/L1Algo/L1Constants.h
+++ b/reco/L1/L1Algo/L1Constants.h
@@ -49,11 +49,6 @@ namespace L1Constants
/// Control flags
namespace control
{
- /// Flag for the radiation length tables usage
- /// true - material budget tables will be used,
- /// false - basic station material info is used
- constexpr bool kIfUseRadLengthTable {true};
-
/// Flag for calling the CAMergeClones procedure ... TODO
constexpr bool kIfMergeClones {true};
diff --git a/reco/L1/L1Algo/L1Fit.h b/reco/L1/L1Algo/L1Fit.h
index e7807c348e..b317ea963b 100644
--- a/reco/L1/L1Algo/L1Fit.h
+++ b/reco/L1/L1Algo/L1Fit.h
@@ -12,7 +12,7 @@
#include "L1Constants.h"
#include "L1Def.h"
-#include "L1MaterialInfo.h"
+#include "L1Material.h"
#include "L1TrackPar.h"
///
@@ -69,13 +69,8 @@ public:
void L1AddMaterial(L1TrackPar& T, fvec radThick, fvec qp0, fvec w);
- void L1AddMaterial(L1TrackPar& T, const L1MaterialInfo& info, fvec qp0, fvec w);
-
-
void L1AddThickMaterial(L1TrackPar& T, fvec radThick, fvec qp0, fvec w, fvec thickness, bool fDownstream);
- void L1AddHalfMaterial(L1TrackPar& T, const L1MaterialInfo& info, fvec qp0);
-
void L1AddPipeMaterial(L1TrackPar& T, fvec qp0, fvec w);
void L1AddTargetMaterial(L1TrackPar& T, fvec qp0, fvec w);
diff --git a/reco/L1/L1Algo/L1FitMaterial.cxx b/reco/L1/L1Algo/L1FitMaterial.cxx
index 2804436807..8491999c77 100644
--- a/reco/L1/L1Algo/L1FitMaterial.cxx
+++ b/reco/L1/L1Algo/L1FitMaterial.cxx
@@ -4,7 +4,7 @@
#include "L1Def.h"
#include "L1Fit.h"
-#include "L1MaterialInfo.h"
+#include "L1Material.h"
#include "L1TrackPar.h"
//#define cnst static const fvec
@@ -375,30 +375,6 @@ void L1Fit::L1AddMaterial(L1TrackPar& T, fvec radThick, fvec qp0, fvec w)
T.C33 += w * (ONE + tyty) * a;
}
-void L1Fit::L1AddMaterial(L1TrackPar& T, const L1MaterialInfo& info, fvec qp0, fvec w)
-{
- cnst ONE = 1.f;
-
- fvec tx = T.tx;
- fvec ty = T.ty;
- fvec txtx = tx * tx;
- fvec tyty = ty * ty;
- fvec txtx1 = txtx + ONE;
- fvec h = txtx + tyty;
- fvec t = sqrt(txtx1 + tyty);
- fvec h2 = h * h;
- fvec qp0t = qp0 * t;
-
- cnst c1 = 0.0136f, c2 = c1 * 0.038f, c3 = c2 * 0.5f, c4 = -c3 / 2.0f, c5 = c3 / 3.0f, c6 = -c3 / 4.0f;
-
- fvec s0 = (c1 + c2 * info.logRadThick + c3 * h + h2 * (c4 + c5 * h + c6 * h2)) * qp0t;
- //fvec a = ( (ONE+fMass2*qp0*qp0t)*info.RadThick*s0*s0 );
- fvec a = ((t + fMass2 * qp0 * qp0t) * info.RadThick * s0 * s0);
-
- T.C22 += w * txtx1 * a;
- T.C32 += w * tx * ty * a;
- T.C33 += w * (ONE + tyty) * a;
-}
// inline void L1Fit::L1AddThickMaterial( L1TrackPar &T, fvec radThick, fvec qp0, fvec thickness=0, fvec w = 1, fvec mass2 = 0.10565f*0.10565f, bool fDownstream = 1 )
// {
@@ -483,31 +459,6 @@ void L1Fit::L1AddThickMaterial(L1TrackPar& T, fvec radThick, fvec qp0, fvec w, f
}
-void L1Fit::L1AddHalfMaterial(L1TrackPar& T, const L1MaterialInfo& info, fvec qp0)
-{
- cnst ONE = 1.f;
- fvec tx = T.tx;
- fvec ty = T.ty;
- fvec txtx = tx * tx;
- fvec tyty = ty * ty;
- fvec txtx1 = txtx + ONE;
- fvec h = txtx + tyty;
- fvec t = sqrt(txtx1 + tyty);
- fvec h2 = h * h;
- fvec qp0t = qp0 * t;
-
- cnst c1(0.0136), c2 = c1 * fvec(0.038), c3 = c2 * fvec(0.5), c4 = -c3 / fvec(2.0), c5 = c3 / fvec(3.0),
- c6 = -c3 / fvec(4.0);
-
- fvec s0 = (c1 + c2 * (info.logRadThick + fvec(log(0.5))) + c3 * h + h2 * (c4 + c5 * h + c6 * h2)) * qp0t;
- //fvec a = ( (ONE+fMass2*qp0*qp0t)*info.RadThick*0.5*s0*s0 );
- fvec a = ((t + fMass2 * qp0 * qp0t) * info.RadThick * fvec(0.5) * s0 * s0);
-
- T.C22 += txtx1 * a;
- T.C32 += tx * ty * a;
- T.C33 += (ONE + tyty) * a;
-}
-
void L1Fit::L1AddPipeMaterial(L1TrackPar& T, fvec qp0, fvec w)
{
cnst ONE = 1.f;
diff --git a/reco/L1/L1Algo/L1InitManager.cxx b/reco/L1/L1Algo/L1InitManager.cxx
index cd35b13995..bd84d8effe 100644
--- a/reco/L1/L1Algo/L1InitManager.cxx
+++ b/reco/L1/L1Algo/L1InitManager.cxx
@@ -42,15 +42,10 @@ void L1InitManager::AddStation(const L1BaseStationInfo& inStation)
L1BaseStationInfo inStationCopy = L1BaseStationInfo(inStation); // make a copy of station so it can be initialized
inStationCopy.SetFieldFunction(fFieldFunction);
- // check, if material map is used
+ // check, if material map is set
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()
- << ". Homogeneous material budget will be used: " << inStationCopy.GetRadThick()[0];
- L1Material material;
- material.SetBins(1, 100);
- material.SetRadThickBin(0, 0, inStationCopy.GetRadThick()[0]);
- inStationCopy.SetMaterialMap(std::move(material));
+ LOG(fatal) << "Station material map was not set for detectorID = "
+ << static_cast<int>(inStationCopy.GetDetectorID()) << ", stationID = " << inStationCopy.GetStationID();
}
// Check station init
diff --git a/reco/L1/L1Algo/L1MaterialInfo.cxx b/reco/L1/L1Algo/L1Material.cxx
similarity index 66%
rename from reco/L1/L1Algo/L1MaterialInfo.cxx
rename to reco/L1/L1Algo/L1Material.cxx
index 442ba6a87d..ada685d20d 100644
--- a/reco/L1/L1Algo/L1MaterialInfo.cxx
+++ b/reco/L1/L1Algo/L1Material.cxx
@@ -2,7 +2,7 @@
SPDX-License-Identifier: GPL-3.0-only
Authors: Sergey Gorbunov, Sergei Zharko [committer] */
-#include "L1MaterialInfo.h"
+#include "L1Material.h"
#include <Logger.h>
@@ -12,60 +12,6 @@
#include "L1Utils.h"
-/************************
- * 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 << "Station thickness (X) [cm]: " << std::setw(12) << std::setfill(' ') << thick[0] << '\n';
- aStream << indent << "Radiation length (X0) [cm]: " << 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();
-}
-
-
-//------------------------------------------------------------------------------------------------------------------------------------
-//
-void L1MaterialInfo::CheckConsistency() const
-{
- /* (i) Checks for the horizontal equality of elements */
- L1Utils::CheckSimdVectorEquality(thick, "L1MaterialInfo::thick");
- L1Utils::CheckSimdVectorEquality(RL, "L1MaterialInfo::RL");
- L1Utils::CheckSimdVectorEquality(RadThick, "L1MaterialInfo::RadThick");
- L1Utils::CheckSimdVectorEquality(logRadThick, "L1MaterialInfo::logRadThick");
-
- /* (ii) Checks for physical sence: thick and RL must be larger then 0. */
- if (thick[0] < fscal(0.)) {
- std::stringstream aStream;
- aStream << "L1MaterialInfo: illegal value for station thickness: (" << thick[0]
- << ", positive value expected) [cm]";
- throw std::logic_error(aStream.str());
- }
-
- if (RL[0] < fscal(0.)) {
- std::stringstream aStream;
- aStream << "L1MaterialInfo: illegal value for station radiation length: (" << RL[0]
- << ", positive value expected) [cm]";
- throw std::logic_error(aStream.str());
- }
-
- /* (iii) Checks for RadThick and logRadThick */
- if (!L1Utils::CmpFloats(RadThick[0] * RL[0], thick[0])) {
- throw std::logic_error(
- "L1MaterialInfo: illegal relation between thickness, radiation length and their ratio (RadThick)");
- }
-
- if (!L1Utils::CmpFloats(std::exp(logRadThick[0]), RadThick[0])) {
- throw std::logic_error("L1MaterialInfo: illegal relation between RadThick and logRadThick data fields");
- }
-}
-
/********************
* L1Material class *
@@ -197,3 +143,79 @@ void L1Material::Swap(L1Material& other) noexcept
std::swap(fFactor, other.fFactor);
std::swap(fTable, other.fTable); // Probably can cause segmentation violation (did not understand)
}
+
+void L1Material::Repare()
+{
+ // correction of the material map: fill empty bins
+ // we assume that bins with radiation thickness == 0. are lacking statistics
+
+ const float cut = 1.e-8;
+ const int n = GetNbins();
+ std::vector<float> oldMap(n * n, 0.);
+
+ bool repeat = 1;
+
+ while (repeat) { // until there are empty bins
+
+ oldMap.assign(oldMap.size(), 0.);
+ for (int iy = 0; iy < n; ++iy) {
+ for (int ix = 0; ix < n; ++ix) {
+ oldMap[iy * n + ix] = GetRadThickBin(ix, iy);
+ }
+ }
+
+ repeat = 0;
+ for (int iy = 0; iy < n; ++iy) {
+ for (int ix = 0; ix < n; ++ix) {
+ if (oldMap[iy * n + ix] >= cut) continue;
+
+ double sum = 0.;
+ double sumw = 0.;
+ // look left
+ for (int i = ix - 1; i >= 0; --i) {
+ float v = oldMap[iy * n + i];
+ if (v >= cut) {
+ double w = 1. / (ix - i);
+ sum += w * v;
+ sumw += w;
+ break;
+ }
+ }
+ // look right
+ for (int i = ix + 1; i < n; ++i) {
+ float v = oldMap[iy * n + i];
+ if (v >= cut) {
+ double w = 1. / (i - ix);
+ sum += w * v;
+ sumw += w;
+ break;
+ }
+ }
+ // look down
+ for (int i = iy - 1; i >= 0; --i) {
+ float v = oldMap[i * n + ix];
+ if (v >= cut) {
+ double w = 1. / (iy - i);
+ sum += w * v;
+ sumw += w;
+ break;
+ }
+ }
+ // look up
+ for (int i = iy + 1; i < n; ++i) {
+ float v = oldMap[i * n + ix];
+ if (v >= cut) {
+ double w = 1. / (i - iy);
+ sum += w * v;
+ sumw += w;
+ break;
+ }
+ }
+ if (sumw > 1.e-8) {
+ SetRadThickBin(ix, iy, sum / sumw);
+ repeat = true;
+ }
+ }
+ }
+ }
+}
diff --git a/reco/L1/L1Algo/L1MaterialInfo.h b/reco/L1/L1Algo/L1Material.h
similarity index 72%
rename from reco/L1/L1Algo/L1MaterialInfo.h
rename to reco/L1/L1Algo/L1Material.h
index 2356f666c1..319be73ebd 100644
--- a/reco/L1/L1Algo/L1MaterialInfo.h
+++ b/reco/L1/L1Algo/L1Material.h
@@ -2,8 +2,8 @@
SPDX-License-Identifier: GPL-3.0-only
Authors: Igor Kulakov, Sergey Gorbunov [committer], Andrey Lebedev, Sergei Zharko */
-#ifndef L1MaterialInfo_h
-#define L1MaterialInfo_h
+#ifndef L1Material_h
+#define L1Material_h
#include <boost/serialization/vector.hpp>
@@ -16,41 +16,6 @@
#include "L1NaN.h"
#include "L1SimdSerializer.h"
-/// 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 {L1NaN::SetNaN<decltype(thick)>()}; ///< Average thickness of the station in arbitary length units
- /// Average radiation length (X0) of the station material in THE SAME UNITS as the thickness
- fvec RL {L1NaN::SetNaN<decltype(RL)>()};
- fvec RadThick {L1NaN::SetNaN<decltype(RadThick)>()}; ///< Average thickness in units of radiation length (X/X0)
- fvec logRadThick {L1NaN::SetNaN<decltype(logRadThick)>()}; ///< Log of average thickness in units of radiation length
-
- /// Verifies class invariant consistency
- void CheckConsistency() const;
-
- /// Checks, if the fields are NaN
- bool IsNaN() const
- {
- return L1NaN::IsNaN(thick) || L1NaN::IsNaN(RL) || L1NaN::IsNaN(RadThick) || L1NaN::IsNaN(logRadThick);
- }
-
- /// String representation of class contents
- /// \param indentLevel number of indent characters in the output
- std::string ToString(int indentLevel = 0) const;
-
- /// Serialization function
- friend class boost::serialization::access;
- template<class Archive>
- void serialize(Archive& ar, const unsigned int)
- {
- ar& thick;
- ar& RL;
- ar& RadThick;
- ar& logRadThick;
- }
-} _fvecalignment;
-
/// Class L1Material describes a map of station thickness in units of radiation length (X0) to the specific point in XY plane
class L1Material {
public:
@@ -116,6 +81,9 @@ public:
/// Swap method
void Swap(L1Material& other) noexcept;
+ /// repare the map - fill empty bins
+ void Repare();
+
private:
int fNbins {L1NaN::SetNaN<decltype(fNbins)>()}; ///< Number of rows (columns) in the material budget table
float fRmax {L1NaN::SetNaN<decltype(fRmax)>()}; ///< Size of the station in x and y dimensions [cm]
diff --git a/reco/L1/L1Algo/L1Parameters.cxx b/reco/L1/L1Algo/L1Parameters.cxx
index 3bd33def76..d2c190ee76 100644
--- a/reco/L1/L1Algo/L1Parameters.cxx
+++ b/reco/L1/L1Algo/L1Parameters.cxx
@@ -134,10 +134,10 @@ void L1Parameters::CheckConsistency() const
for (int iSt = 0; iSt < fNstationsActiveTotal; ++iSt) {
fStations[iSt].CheckConsistency();
- if (fStations[iSt].z[0] < fTargetPos[2][0]) {
+ 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].z[0] << " [cm] < z_targ = " << fTargetPos[2][0] << " [cm])";
+ << "(z_st = " << fStations[iSt].fZ[0] << " [cm] < z_targ = " << fTargetPos[2][0] << " [cm])";
throw std::logic_error(msg.str());
}
}
diff --git a/reco/L1/L1Algo/L1Parameters.h b/reco/L1/L1Algo/L1Parameters.h
index cd78ec2104..4fa9318884 100644
--- a/reco/L1/L1Algo/L1Parameters.h
+++ b/reco/L1/L1Algo/L1Parameters.h
@@ -17,7 +17,7 @@
#include "L1CAIteration.h"
#include "L1Constants.h"
-#include "L1MaterialInfo.h"
+#include "L1Material.h"
#include "L1Station.h"
#include "L1Utils.h"
#include "L1Vector.h"
diff --git a/reco/L1/L1Algo/L1Station.cxx b/reco/L1/L1Algo/L1Station.cxx
index 8c3bf9905b..d8f46fa761 100644
--- a/reco/L1/L1Algo/L1Station.cxx
+++ b/reco/L1/L1Algo/L1Station.cxx
@@ -41,7 +41,7 @@ void L1Station::CheckConsistency() const
* SIMD vector checks: all the words in a SIMD vector must be equal
*/
- L1Utils::CheckSimdVectorEquality(z, "L1Station::z");
+ L1Utils::CheckSimdVectorEquality(fZ, "L1Station::fZ");
L1Utils::CheckSimdVectorEquality(Rmin, "L1Station::Rmin");
L1Utils::CheckSimdVectorEquality(Rmax, "L1Station::Rmax");
L1Utils::CheckSimdVectorEquality(dt, "L1Station::dt");
@@ -113,16 +113,16 @@ std::string L1Station::ToString(int verbosityLevel, int indentLevel) const
std::stringstream aStream {};
constexpr char indentChar = '\t';
std::string indent(indentLevel, indentChar);
- if (verbosityLevel == 0) { aStream << "station type = " << type << ", z = " << z[0] << " cm"; }
+ if (verbosityLevel == 0) { aStream << "station type = " << type << ", z = " << fZ[0] << " cm"; }
else {
aStream << '\n';
- aStream << indent << "z pos [cm]: " << std::setw(12) << std::setfill(' ') << z[0] << '\n';
+ aStream << indent << "z pos [cm]: " << std::setw(12) << std::setfill(' ') << fZ[0] << '\n';
+ aStream << indent << "z thick [cm]: " << std::setw(12) << std::setfill(' ') << fZthick[0] << '\n';
aStream << indent << "R_min [cm]: " << std::setw(12) << std::setfill(' ') << Rmin[0] << '\n';
aStream << indent << "R_max [cm]: " << std::setw(12) << std::setfill(' ') << Rmax[0] << '\n';
aStream << indent << "Station type: " << std::setw(12) << std::setfill(' ') << type << '\n';
aStream << indent << "Is time used: " << std::setw(12) << std::setfill(' ') << timeInfo << '\n';
aStream << indent << "Is in field: " << std::setw(12) << std::setfill(' ') << fieldStatus << '\n';
- aStream << materialInfo.ToString(indentLevel) << '\n';
if (verbosityLevel > 3) {
aStream << indent << "Field approcimation coefficients:\n";
diff --git a/reco/L1/L1Algo/L1Station.h b/reco/L1/L1Algo/L1Station.h
index c0103b1b7f..035d99974c 100644
--- a/reco/L1/L1Algo/L1Station.h
+++ b/reco/L1/L1Algo/L1Station.h
@@ -9,7 +9,7 @@
#include <type_traits>
#include "L1Field.h"
-#include "L1MaterialInfo.h"
+#include "L1Material.h"
#include "L1NaN.h"
#include "L1SimdSerializer.h"
#include "L1UMeasurementInfo.h"
@@ -25,12 +25,11 @@ public:
int timeInfo {L1NaN::SetNaN<decltype(timeInfo)>()}; ///< flag: if time information can be used
int fieldStatus {L1NaN::SetNaN<decltype(
fieldStatus)>()}; ///< flag: 1 - station is INSIDE the field, 0 - station is OUTSIDE the field (replace with enum)
- fvec z {L1NaN::SetNaN<decltype(z)>()}; ///< z position of station [cm]
- fvec Rmin {L1NaN::SetNaN<decltype(Rmin)>()}; ///< min radius of the station [cm]
- fvec Rmax {L1NaN::SetNaN<decltype(Rmax)>()}; ///< max radius of the station [cm]
- fvec dt {L1NaN::SetNaN<decltype(dt)>()}; ///< time resolution [ns]
- /// structure containing station thickness(X), rad. length (X0), X/X0 and log(X/X0) values
- L1MaterialInfo materialInfo {};
+ fvec fZ {L1NaN::SetNaN<decltype(fZ)>()}; ///< z position of station [cm]
+ fvec fZthick {L1NaN::SetNaN<decltype(fZthick)>()}; ///< z thickness of the station [cm]
+ fvec Rmin {L1NaN::SetNaN<decltype(Rmin)>()}; ///< min radius of the station [cm]
+ fvec Rmax {L1NaN::SetNaN<decltype(Rmax)>()}; ///< max radius of the station [cm]
+ fvec dt {L1NaN::SetNaN<decltype(dt)>()}; ///< time resolution [ns]
L1FieldSlice fieldSlice {};
L1UMeasurementInfo frontInfo {};
L1UMeasurementInfo backInfo {};
@@ -47,12 +46,12 @@ public:
ar& timeInfo;
ar& fieldStatus;
- ar& z;
+ ar& fZ;
+ ar& fZthick;
ar& Rmin;
ar& Rmax;
ar& dt;
- ar& materialInfo;
ar& fieldSlice;
ar& frontInfo;
ar& backInfo;
diff --git a/reco/L1/L1Algo/L1TrackFitter.cxx b/reco/L1/L1Algo/L1TrackFitter.cxx
index 60c1b48d72..79028053cb 100644
--- a/reco/L1/L1Algo/L1TrackFitter.cxx
+++ b/reco/L1/L1Algo/L1TrackFitter.cxx
@@ -109,9 +109,9 @@ void L1Algo::KFTrackFitter_simple() // TODO: Add pipe.
// static L1FieldRegion fld _fvecalignment;
L1FieldValue fldB0, fldB1, fldB2 _fvecalignment;
L1FieldRegion fld _fvecalignment;
- fvec fldZ0 = sta1.z; // suppose field is smoth
- fvec fldZ1 = sta2.z;
- fvec fldZ2 = sta0.z;
+ fvec fldZ0 = sta1.fZ; // suppose field is smoth
+ fvec fldZ1 = sta2.fZ;
+ fvec fldZ2 = sta0.fZ;
sta1.fieldSlice.GetFieldValue(x1, y1, fldB0);
sta2.fieldSlice.GetFieldValue(x2, y2, fldB1);
@@ -132,12 +132,7 @@ void L1Algo::KFTrackFitter_simple() // TODO: Add pipe.
L1ExtrapolateLine(T, hit.z);
// T.L1Extrapolate( sta.z, qp0, fld );
// L1Extrapolate( T, hit.z, qp0, fld );
- if constexpr (L1Constants::control::kIfUseRadLengthTable) {
- fit.L1AddMaterial(T, fParameters.GetMaterialThickness(i, T.x, T.y), qp0, ONE);
- }
- else {
- fit.L1AddMaterial(T, sta.materialInfo, qp0, ONE);
- }
+ fit.L1AddMaterial(T, fParameters.GetMaterialThickness(i, T.x, T.y), qp0, ONE);
// if (ista == fNstationsBeforePipe - 1) fit.L1AddPipeMaterial( T, qp0, 1);
@@ -152,7 +147,7 @@ void L1Algo::KFTrackFitter_simple() // TODO: Add pipe.
fldZ1 = fldZ2;
sta.fieldSlice.GetFieldValue(x, y, fldB2);
- fldZ2 = sta.z;
+ fldZ2 = sta.fZ;
fld.Set(fldB2, fldZ2, fldB1, fldZ1, fldB0, fldZ0);
} // i
@@ -242,9 +237,9 @@ void L1Algo::KFTrackFitter_simple() // TODO: Add pipe.
// static L1FieldRegion fld _fvecalignment;
L1FieldValue fldB0, fldB1, fldB2 _fvecalignment;
L1FieldRegion fld _fvecalignment;
- fvec fldZ0 = sta1.z;
- fvec fldZ1 = sta2.z;
- fvec fldZ2 = sta0.z;
+ fvec fldZ0 = sta1.fZ;
+ fvec fldZ1 = sta2.fZ;
+ fvec fldZ2 = sta0.fZ;
sta1.fieldSlice.GetFieldValue(x1, y1, fldB0);
sta2.fieldSlice.GetFieldValue(x2, y2, fldB1);
@@ -265,12 +260,8 @@ void L1Algo::KFTrackFitter_simple() // TODO: Add pipe.
L1ExtrapolateLine(T, hit.z);
// T.L1Extrapolate( sta.z, qp0, fld );
// L1Extrapolate( T, hit.z, qp0, fld );
- if constexpr (L1Constants::control::kIfUseRadLengthTable) {
- fit.L1AddMaterial(T, fParameters.GetMaterialThickness(i, T.x, T.y), qp0, ONE);
- }
- else {
- fit.L1AddMaterial(T, sta.materialInfo, qp0, ONE);
- }
+ fit.L1AddMaterial(T, fParameters.GetMaterialThickness(i, T.x, T.y), qp0, ONE);
+
// if (ista == fNstationsBeforePipe) fit.L1AddPipeMaterial( T, qp0, 1);
L1Filter(T, sta.frontInfo, u);
L1Filter(T, sta.backInfo, v);
@@ -280,7 +271,7 @@ void L1Algo::KFTrackFitter_simple() // TODO: Add pipe.
fldZ0 = fldZ1;
fldZ1 = fldZ2;
sta.fieldSlice.GetFieldValue(x, y, fldB2);
- fldZ2 = sta.z;
+ fldZ2 = sta.fZ;
fld.Set(fldB2, fldZ2, fldB1, fldZ1, fldB0, fldZ0);
}
@@ -396,7 +387,7 @@ void L1Algo::L1KFTrackFitter()
fvec ZSta[L1Constants::size::kMaxNstations];
for (int ista = 0; ista < nStations; ista++) {
- ZSta[ista] = sta[ista].z;
+ ZSta[ista] = sta[ista].fZ;
}
unsigned short N_vTracks = fTracks.size(); // number of tracks processed per one SSE register
@@ -421,6 +412,8 @@ void L1Algo::L1KFTrackFitter()
z[ista] = ZSta[ista];
}
+ //fmask isFieldPresent = fmask::Zero();
+
for (int iVec = 0; iVec < nTracks_SIMD; iVec++) {
int nHitsTrack = t[iVec]->NHits;
@@ -430,8 +423,11 @@ void L1Algo::L1KFTrackFitter()
const L1Hit& hit = fInputData.GetHit(fRecoHits[start_hit++]);
const int ista = hit.iSt;
- iSta[ih] = ista;
- w[ista][iVec] = 1.;
+
+ //if (sta[ista].fieldStatus) { isFieldPresent[iVec] = true; }
+
+ iSta[ih] = ista;
+ w[ista][iVec] = 1.;
if (sta[ista].timeInfo) { w_time[ista][iVec] = 1.; }
u[ista][iVec] = hit.u;
@@ -497,6 +493,8 @@ void L1Algo::L1KFTrackFitter()
if (kGlobal == fTrackingMode || kMcbm == fTrackingMode) { tr.qp = fvec(0.); }
+ //tr.qp = iif(isFieldPresent, tr.qp, fvec(1. / 0.25));
+
for (int iter = 0; iter < 2; iter++) { // 1.5 iterations
fvec qp01 = tr.qp;
@@ -543,13 +541,8 @@ void L1Algo::L1KFTrackFitter()
//fit.AddPipeMaterial(qp01, wExtr);
//fit.EnergyLossCorrection(fit.fPipeRadThick, qp01, fvec(1.f), wExtr);
//}
- if constexpr (L1Constants::control::kIfUseRadLengthTable) {
- fit.AddMaterial(fParameters.GetMaterialThickness(ista, tr.x, tr.y), qp01, wExtr);
- fit.EnergyLossCorrection(fParameters.GetMaterialThickness(ista, tr.x, tr.y), qp01, fvec(1.f), wExtr);
- }
- else {
- fit.AddMaterial(sta[ista].materialInfo, qp01, wExtr);
- }
+ fit.AddMaterial(fParameters.GetMaterialThickness(ista, tr.x, tr.y), qp01, wExtr);
+ fit.EnergyLossCorrection(fParameters.GetMaterialThickness(ista, tr.x, tr.y), qp01, fvec(1.f), wExtr);
fit.Filter(sta[ista].frontInfo, u[ista], du2[ista], w1);
fit.Filter(sta[ista].backInfo, v[ista], dv2[ista], w1);
@@ -700,13 +693,8 @@ void L1Algo::L1KFTrackFitter()
//fit.AddPipeMaterial(qp01, wExtr);
//fit.EnergyLossCorrection(fit.fPipeRadThick, qp01, fvec(-1.f), wExtr);
//}
- if constexpr (L1Constants::control::kIfUseRadLengthTable) {
- fit.AddMaterial(fParameters.GetMaterialThickness(ista, tr.x, tr.y), qp01, wExtr);
- fit.EnergyLossCorrection(fParameters.GetMaterialThickness(ista, tr.x, tr.y), qp01, fvec(-1.f), wExtr);
- }
- else {
- fit.AddMaterial(sta[ista].materialInfo, qp01, wExtr);
- }
+ fit.AddMaterial(fParameters.GetMaterialThickness(ista, tr.x, tr.y), qp01, wExtr);
+ fit.EnergyLossCorrection(fParameters.GetMaterialThickness(ista, tr.x, tr.y), qp01, fvec(-1.f), wExtr);
fit.Filter(sta[ista].frontInfo, u[ista], du2[ista], w1);
fit.Filter(sta[ista].backInfo, v[ista], dv2[ista], w1);
diff --git a/reco/L1/L1Algo/L1TrackParFit.cxx b/reco/L1/L1Algo/L1TrackParFit.cxx
index 974c6e10ff..d19f8c0149 100644
--- a/reco/L1/L1Algo/L1TrackParFit.cxx
+++ b/reco/L1/L1Algo/L1TrackParFit.cxx
@@ -885,33 +885,6 @@ void L1TrackParFit::AddThickMaterial(fvec radThick, fvec qp0, fvec w, fvec thick
}
-void L1TrackParFit::AddMaterial(const L1MaterialInfo& info, fvec qp0, fvec w)
-{
- cnst ONE = 1.f;
-
- fvec tx = fTr.tx;
- fvec ty = fTr.ty;
- // fvec time = fTr.t;
- fvec txtx = tx * tx;
- fvec tyty = ty * ty;
- fvec txtx1 = txtx + ONE;
- fvec h = txtx + tyty;
- fvec t = sqrt(txtx1 + tyty);
- fvec h2 = h * h;
- fvec qp0t = qp0 * t;
-
- cnst c1 = 0.0136f, c2 = c1 * 0.038f, c3 = c2 * 0.5f, c4 = -c3 / 2.0f, c5 = c3 / 3.0f, c6 = -c3 / 4.0f;
-
- fvec s0 = (c1 + c2 * info.logRadThick + c3 * h + h2 * (c4 + c5 * h + c6 * h2)) * qp0t;
- //fvec a = ( (ONE+mass2*qp0*qp0t)*info.RadThick*s0*s0 );
- fvec a = ((t + fMass2 * qp0 * qp0t) * info.RadThick * s0 * s0);
-
- fTr.C22 += w * txtx1 * a;
- fTr.C32 += w * tx * ty * a;
- fTr.C33 += w * (ONE + tyty) * a;
-}
-
-
void L1TrackParFit::EnergyLossCorrection(const fvec& radThick, fvec& qp0, fvec direction, fvec w)
{
const fvec qp2cut(1. / (10. * 10.)); // 10 GeV cut
diff --git a/reco/L1/L1Algo/L1TrackParFit.h b/reco/L1/L1Algo/L1TrackParFit.h
index 2a800cf265..20caaa3290 100644
--- a/reco/L1/L1Algo/L1TrackParFit.h
+++ b/reco/L1/L1Algo/L1TrackParFit.h
@@ -7,7 +7,6 @@
#include "L1Def.h"
#include "L1Field.h"
-#include "L1MaterialInfo.h"
#include "L1TrackPar.h"
#include "L1UMeasurementInfo.h"
#include "L1XYMeasurementInfo.h"
@@ -70,8 +69,6 @@ public:
void EnergyLossCorrectionCarbon(const fvec& radThick, fvec& qp0, fvec direction, fvec w);
void EnergyLossCorrectionAl(const fvec& radThick, fvec& qp0, fvec direction, fvec w);
- void AddMaterial(const L1MaterialInfo& info, fvec qp0, fvec w);
-
void AddMaterial(const fvec& radThick, fvec qp0, fvec w = 1);
void AddThickMaterial(fvec radThick, fvec qp0, fvec w, fvec thickness, bool fDownstream);
diff --git a/reco/L1/L1Algo/utils/L1AlgoDraw.cxx b/reco/L1/L1Algo/utils/L1AlgoDraw.cxx
index 43faddf73f..e03940db1b 100644
--- a/reco/L1/L1Algo/utils/L1AlgoDraw.cxx
+++ b/reco/L1/L1Algo/utils/L1AlgoDraw.cxx
@@ -593,7 +593,7 @@ void L1AlgoDraw::DrawInputHits()
TLine* line = new TLine();
line->SetLineColor(StaColor);
- line->DrawLine(st.z[0], -st.Rmax[0], st.z[0], st.Rmax[0]);
+ line->DrawLine(st.fZ[0], -st.Rmax[0], st.fZ[0], st.Rmax[0]);
TPolyMarker* pmyz = new TPolyMarker(n_poly, z_poly, y_poly);
pmyz->SetMarkerColor(mcolor[ista]);
@@ -609,7 +609,7 @@ void L1AlgoDraw::DrawInputHits()
XZ->cd();
- line->DrawLine(st.z[0], -st.Rmax[0], st.z[0], st.Rmax[0]);
+ line->DrawLine(st.fZ[0], -st.Rmax[0], st.fZ[0], st.Rmax[0]);
TPolyMarker* pmxz = new TPolyMarker(n_poly, z_poly, x_poly);
pmxz->SetMarkerColor(mcolor[ista]);
@@ -714,7 +714,7 @@ void L1AlgoDraw::DrawRestHits(L1HitIndex_t* StsRestHitsStartIndex, L1HitIndex_t*
TLine* line = new TLine();
line->SetLineColor(StaColor);
- line->DrawLine(st.z[0], -st.Rmax[0], st.z[0], st.Rmax[0]);
+ line->DrawLine(st.fZ[0], -st.Rmax[0], st.fZ[0], st.Rmax[0]);
TPolyMarker* pmyz = new TPolyMarker(n_poly, z_poly, y_poly);
pmyz->SetMarkerColor(mcolor[ista]);
@@ -730,7 +730,7 @@ void L1AlgoDraw::DrawRestHits(L1HitIndex_t* StsRestHitsStartIndex, L1HitIndex_t*
XZ->cd();
- line->DrawLine(st.z[0], -st.Rmax[0], st.z[0], st.Rmax[0]);
+ line->DrawLine(st.fZ[0], -st.Rmax[0], st.fZ[0], st.Rmax[0]);
TPolyMarker* pmxz = new TPolyMarker(n_poly, z_poly, x_poly);
pmxz->SetMarkerColor(mcolor[ista]);
diff --git a/reco/L1/ParticleFinder/CbmL1PFFitter.cxx b/reco/L1/ParticleFinder/CbmL1PFFitter.cxx
index d2aa8055a8..c55858e65d 100644
--- a/reco/L1/ParticleFinder/CbmL1PFFitter.cxx
+++ b/reco/L1/ParticleFinder/CbmL1PFFitter.cxx
@@ -42,7 +42,7 @@
#include "L1Field.h"
#include "L1Filtration.h"
#include "L1Fit.h"
-#include "L1MaterialInfo.h"
+#include "L1Material.h"
#include "L1Station.h"
#include "L1TrackPar.h"
@@ -206,7 +206,7 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo)
// get hits of current track
for (i = 0; i < nHits; i++) {
w[i] = ZERO;
- z[i] = sta[i].z;
+ z[i] = sta[i].fZ;
}
for (iVec = 0; iVec < nTracks_SIMD; iVec++) {
@@ -300,14 +300,9 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo)
fit.L1AddPipeMaterial(T, qp0, wIn);
fit.EnergyLossCorrection(T, fit.PipeRadThick, qp0, fvec(-1.f), wIn);
}
- if constexpr (L1Constants::control::kIfUseRadLengthTable) {
- fit.L1AddMaterial(T, CbmL1::Instance()->fpAlgo->GetParameters()->GetMaterialThickness(i, T.x, T.y), qp0, wIn);
- fit.EnergyLossCorrection(T, CbmL1::Instance()->fpAlgo->GetParameters()->GetMaterialThickness(i, T.x, T.y), qp0,
- -1, wIn);
- }
- else {
- fit.L1AddMaterial(T, sta[i].materialInfo, qp0, wIn);
- }
+ fit.L1AddMaterial(T, CbmL1::Instance()->fpAlgo->GetParameters()->GetMaterialThickness(i, T.x, T.y), qp0, wIn);
+ fit.EnergyLossCorrection(T, CbmL1::Instance()->fpAlgo->GetParameters()->GetMaterialThickness(i, T.x, T.y), qp0,
+ -1, wIn);
L1Filter(T, sta[i].frontInfo, u[i], w1);
L1Filter(T, sta[i].backInfo, v[i], w1);
@@ -377,14 +372,9 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo)
fit.L1AddPipeMaterial(T, qp0, wIn);
fit.EnergyLossCorrection(T, fit.PipeRadThick, qp0, fvec(1.f), wIn);
}
- if constexpr (L1Constants::control::kIfUseRadLengthTable) {
- fit.L1AddMaterial(T, CbmL1::Instance()->fpAlgo->GetParameters()->GetMaterialThickness(i, T.x, T.y), qp0, wIn);
- fit.EnergyLossCorrection(T, CbmL1::Instance()->fpAlgo->GetParameters()->GetMaterialThickness(i, T.x, T.y), qp0,
- 1, wIn);
- }
- else {
- fit.L1AddMaterial(T, sta[i].materialInfo, qp0, wIn);
- }
+ fit.L1AddMaterial(T, CbmL1::Instance()->fpAlgo->GetParameters()->GetMaterialThickness(i, T.x, T.y), qp0, wIn);
+ fit.EnergyLossCorrection(T, CbmL1::Instance()->fpAlgo->GetParameters()->GetMaterialThickness(i, T.x, T.y), qp0, 1,
+ wIn);
L1Filter(T, sta[i].frontInfo, u[i], w1);
L1Filter(T, sta[i].backInfo, v[i], w1);
@@ -450,7 +440,7 @@ void CbmL1PFFitter::GetChiToVertex(vector<CbmStsTrack>& Tracks, vector<PFFieldRe
const L1Station* sta = CbmL1::Instance()->fpAlgo->GetParameters()->GetStations().begin();
fvec* zSta = new fvec[nStations];
for (int iSta = 0; iSta < nStations; iSta++) {
- zSta[iSta] = sta[iSta].z;
+ zSta[iSta] = sta[iSta].fZ;
}
field.reserve(Tracks.size());
--
GitLab