From 890f0b2a16ee5c374a7cd633258584dc885adb10 Mon Sep 17 00:00:00 2001
From: "se.gorbunov" <se.gorbunov@gsi.de>
Date: Tue, 11 Oct 2022 21:35:00 +0000
Subject: [PATCH] L1: dev option to use the original field
---
reco/L1/CbmL1.cxx | 3 +-
reco/L1/CbmL1Performance.cxx | 82 +++---------------------------
reco/L1/L1Algo/L1Algo.cxx | 2 +
reco/L1/L1Algo/L1Extrapolation.cxx | 18 +++----
reco/L1/L1Algo/L1Field.cxx | 40 ++-------------
reco/L1/L1Algo/L1Field.h | 22 ++++----
reco/L1/L1Algo/L1InitManager.h | 3 ++
reco/L1/L1Algo/L1Parameters.h | 4 ++
8 files changed, 41 insertions(+), 133 deletions(-)
diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx
index 492ad9fba2..92bad00f2a 100644
--- a/reco/L1/CbmL1.cxx
+++ b/reco/L1/CbmL1.cxx
@@ -215,7 +215,6 @@ InitStatus CbmL1::Init()
fInitManager.DevSetIgnoreHitSearchAreas(true);
}
-
if (L1Algo::TrackingMode::kGlobal == fTrackingMode) {
//at the moment trd2d tracking only
fUseMVD = 0;
@@ -224,7 +223,7 @@ InitStatus CbmL1::Init()
fUseTRD = 1;
fUseTOF = 0;
- L1FieldRegion::gkUseOriginalField = true;
+ fInitManager.DevSetUseOfOriginalField();
fInitManager.DevSetIgnoreHitSearchAreas(true);
//fInitManager.DevSetFitSingletsFromTarget(true);
//fInitManager.DevSetIsMatchDoubletsViaMc(true);
diff --git a/reco/L1/CbmL1Performance.cxx b/reco/L1/CbmL1Performance.cxx
index d073ee2d09..189ab0e573 100644
--- a/reco/L1/CbmL1Performance.cxx
+++ b/reco/L1/CbmL1Performance.cxx
@@ -1229,25 +1229,10 @@ void CbmL1::TrackFitPerformance()
L1TrackPar trPar(it->T, it->C);
L1FieldRegion fld _fvecalignment;
- L1FieldValue B[3], targB _fvecalignment;
- float z[3] = {0.f, 0.f, 0.f};
- int ih = 0;
- for (unsigned int iMCPoint = 0; iMCPoint < mc.Points.size(); iMCPoint++) {
- const int iMCP = mc.Points[iMCPoint];
- CbmL1MCPoint& mcP = fvMCPoints[iMCP];
- const L1Station& st = fpAlgo->GetParameters()->GetStation(mcP.iStation);
- z[ih] = st.z[0];
- if (ih > 0 && (z[ih] - z[ih - 1]) < 0.1) continue;
- st.fieldSlice.GetFieldValue(mcP.x, mcP.y, B[ih]);
- ih++;
- if (ih == 3) break;
- }
- if (ih < 3) continue;
+ fld.SetUseOriginalField();
CbmL1MCPoint& mcP = fvMCPoints[mc.Points[0]];
- targB = fpAlgo->GetParameters()->GetVertexFieldValue();
- fld.Set(B[0], z[0], B[1], z[1], B[2], z[2]);
L1Extrapolate(trPar, mcP.zIn, trPar.qp, fld);
const L1TrackPar& tr = trPar;
@@ -1345,23 +1330,9 @@ void CbmL1::TrackFitPerformance()
CbmL1MCTrack mc = *(it->GetMCTracks()[0]);
L1TrackPar trPar(it->TLast, it->CLast);
L1FieldRegion fld _fvecalignment;
- L1FieldValue B[3], targB _fvecalignment;
- float z[3] = {0.f, 0.f, 0.f};
- int ih = 0;
- for (unsigned int iMCPoint = 0; iMCPoint < mc.Points.size(); iMCPoint++) {
- const int iMCP = mc.Points[iMCPoint];
- CbmL1MCPoint& mcP = fvMCPoints[iMCP];
- const L1Station& st = fpAlgo->GetParameters()->GetStation(mcP.iStation);
- z[ih] = st.z[0];
- if (ih > 0 && (z[ih] - z[ih - 1]) < 0.1) continue;
- st.fieldSlice.GetFieldValue(mcP.x, mcP.y, B[ih]);
- ih++;
- if (ih == 3) break;
- }
- if (ih < 3) continue;
+ fld.SetUseOriginalField();
+
CbmL1MCPoint& mcP = fvMCPoints[iMC];
- targB = fpAlgo->GetParameters()->GetVertexFieldValue();
- fld.Set(B[0], z[0], B[1], z[1], B[2], z[2]);
L1Extrapolate(trPar, mcP.zOut, trPar.qp, fld);
const L1TrackPar& tr = trPar;
@@ -1419,18 +1390,7 @@ void CbmL1::TrackFitPerformance()
{ // extrapolate to vertex
L1FieldRegion fld _fvecalignment;
- L1FieldValue B[3] _fvecalignment;
- float z[3] = {0.f, 0.f, 0.f};
- for (unsigned int ih = 0; ih < 3; ih++) {
- if (ih >= mc.Points.size()) continue; //If nofMCPoints in track < 3
- const int iMCP = mc.Points[ih];
- CbmL1MCPoint& mcP = fvMCPoints[iMCP];
- const L1Station& st = fpAlgo->GetParameters()->GetStation(mcP.iStation);
- z[ih] = st.z[0];
- st.fieldSlice.GetFieldValue(mcP.x, mcP.y, B[ih]);
- };
- fld.Set(B[0], z[0], B[1], z[1], B[2], z[2]);
-
+ fld.SetUseOriginalField();
L1Extrapolate(trPar, mc.z, trPar.qp, fld);
// add material
const int fSta = fvHitStore[it->Hits[0]].iStation;
@@ -1465,7 +1425,7 @@ void CbmL1::TrackFitPerformance()
h_fitSV[1]->Fill((tr.y[0] - mc.y));
h_fitSV[2]->Fill((tr.tx[0] - mc.px / mc.pz) * 1.e3);
h_fitSV[3]->Fill((tr.ty[0] - mc.py / mc.pz) * 1.e3);
- h_fitSV[4]->Fill(fabs(1. / tr.qp[0]) / mc.p - 1);
+ h_fitSV[4]->Fill(fabs(1. / tr.qp[0]) / mc.p - 1.);
if (std::isfinite(tr.C00[0]) && tr.C00[0] > 0) h_fitSV[5]->Fill((tr.x[0] - mc.x) / sqrt(tr.C00[0]));
if (std::isfinite(tr.C11[0]) && tr.C11[0] > 0) h_fitSV[6]->Fill((tr.y[0] - mc.y) / sqrt(tr.C11[0]));
if (std::isfinite(tr.C22[0]) && tr.C22[0] > 0) h_fitSV[7]->Fill((tr.tx[0] - mc.px / mc.pz) / sqrt(tr.C22[0]));
@@ -1482,21 +1442,7 @@ void CbmL1::TrackFitPerformance()
#ifdef L1EXTRAPOLATE
{ // extrapolate to vertex
L1FieldRegion fld _fvecalignment;
- L1FieldValue B[3], targB _fvecalignment;
- float z[3];
-
- targB = fpAlgo->GetParameters()->GetVertexFieldValue();
-
- int ih = 1;
- for (unsigned int iHit = 0; iHit < it->Hits.size(); iHit++) {
- const int iStation = fvHitStore[it->Hits[iHit]].iStation;
- const L1Station& st = fpAlgo->GetParameters()->GetStation(iStation);
- z[ih] = st.z[0];
- st.fieldSlice.GetFieldValue(fvHitStore[it->Hits[iHit]].x, fvHitStore[it->Hits[iHit]].y, B[ih]);
- ih++;
- if (ih == 3) break;
- }
- if (ih < 3) continue;
+ fld.SetUseOriginalField();
// add material
const int fSta = fvHitStore[it->Hits[0]].iStation;
@@ -1509,30 +1455,16 @@ void CbmL1::TrackFitPerformance()
&& (dir * (mc.z - fpAlgo->GetParameters()->GetStation(iSta).z[0]) > 0);
iSta += dir) {
- z[0] = fpAlgo->GetParameters()->GetStation(iSta).z[0];
- float dz = z[1] - z[0];
- fpAlgo->GetParameters()->GetStation(iSta).fieldSlice.GetFieldValue(trPar.x - trPar.tx * dz,
- trPar.y - trPar.ty * dz, B[0]);
- fld.Set(B[0], z[0], B[1], z[1], B[2], z[2]);
-
L1Extrapolate(trPar, fpAlgo->GetParameters()->GetStation(iSta).z[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),
- trPar.qp, fvec(1.f), fvec(1.f));
+ trPar.qp, fvec(1.), fvec(1.));
if (iSta + dir == fNMvdStations - 1) {
fit.L1AddPipeMaterial(trPar, trPar.qp, 1);
fit.EnergyLossCorrection(trPar, fit.PipeRadThick, trPar.qp, fvec(1.f), fvec(1.f));
}
- B[2] = B[1];
- z[2] = z[1];
- B[1] = B[0];
- z[1] = z[0];
}
-
- z[0] = fpAlgo->GetParameters()->GetTargetPositionZ()[0];
- B[0] = targB;
- fld.Set(B[0], z[0], B[1], z[1], B[2], z[2]);
L1Extrapolate(trPar, mc.z, trPar.qp, fld);
}
if (mc.z != trPar.z[0]) continue;
diff --git a/reco/L1/L1Algo/L1Algo.cxx b/reco/L1/L1Algo/L1Algo.cxx
index 2b0375f694..065734e07e 100644
--- a/reco/L1/L1Algo/L1Algo.cxx
+++ b/reco/L1/L1Algo/L1Algo.cxx
@@ -149,6 +149,8 @@ void L1Algo::ReceiveParameters(L1Parameters&& parameters)
fGhostSuppression = fParameters.GetGhostSuppression();
fMomentumCutOff = fParameters.GetMomentumCutOff();
+ L1FieldRegion::ForceUseOfOriginalField(fParameters.DevIsUseOfOriginalField());
+
LOG(info) << fParameters.ToString(3);
}
diff --git a/reco/L1/L1Algo/L1Extrapolation.cxx b/reco/L1/L1Algo/L1Extrapolation.cxx
index cbb081f569..99ea46080f 100644
--- a/reco/L1/L1Algo/L1Extrapolation.cxx
+++ b/reco/L1/L1Algo/L1Extrapolation.cxx
@@ -352,18 +352,16 @@ void L1ExtrapolateRungeKutta // extrapolates track parameters and returns jacob
int step;
int i;
- fvec B[4][3];
for (step = 0; step < 4; ++step) {
- F.Get(z_in + a[step] * h, B[step]);
- }
- for (step = 0; step < 4; ++step) {
for (i = 0; i < 4; ++i) {
if (step == 0) { x[i] = x0[i]; }
else {
x[i] = x0[i] + b[step] * k[step * 4 - 4 + i];
}
}
+ fvec B[3];
+ F.Get(x[0], x[1], z_in + a[step] * h, B);
fvec tx = x[2];
fvec ty = x[3];
@@ -377,13 +375,13 @@ void L1ExtrapolateRungeKutta // extrapolates track parameters and returns jacob
// fvec I_tx2ty2 = qp0 * hC / tx2ty2 ; unsused ???
tx2ty2 *= hC;
fvec tx2ty2qp = tx2ty2 * qp0;
- Ax[step] = (txty * B[step][0] + ty * B[step][2] - (fvec(1.) + tx2) * B[step][1]) * tx2ty2;
- Ay[step] = (-txty * B[step][1] - tx * B[step][2] + (fvec(1.) + ty2) * B[step][0]) * tx2ty2;
+ Ax[step] = (txty * B[0] + ty * B[2] - (fvec(1.) + tx2) * B[1]) * tx2ty2;
+ Ay[step] = (-txty * B[1] - tx * B[2] + (fvec(1.) + ty2) * B[0]) * tx2ty2;
- Ax_tx[step] = Ax[step] * tx * I_tx2ty21 + (ty * B[step][0] - fvec(2.) * tx * B[step][1]) * tx2ty2qp;
- Ax_ty[step] = Ax[step] * ty * I_tx2ty21 + (tx * B[step][0] + B[step][2]) * tx2ty2qp;
- Ay_tx[step] = Ay[step] * tx * I_tx2ty21 + (-ty * B[step][1] - B[step][2]) * tx2ty2qp;
- Ay_ty[step] = Ay[step] * ty * I_tx2ty21 + (-tx * B[step][1] + fvec(2.) * ty * B[step][0]) * tx2ty2qp;
+ Ax_tx[step] = Ax[step] * tx * I_tx2ty21 + (ty * B[0] - fvec(2.) * tx * B[1]) * tx2ty2qp;
+ Ax_ty[step] = Ax[step] * ty * I_tx2ty21 + (tx * B[0] + B[2]) * tx2ty2qp;
+ Ay_tx[step] = Ay[step] * tx * I_tx2ty21 + (-ty * B[1] - B[2]) * tx2ty2qp;
+ Ay_ty[step] = Ay[step] * ty * I_tx2ty21 + (-tx * B[1] + fvec(2.) * ty * B[0]) * tx2ty2qp;
step4 = step * 4;
k[step4] = tx * h;
diff --git a/reco/L1/L1Algo/L1Field.cxx b/reco/L1/L1Algo/L1Field.cxx
index 9f4091592b..c91cf325b3 100644
--- a/reco/L1/L1Algo/L1Field.cxx
+++ b/reco/L1/L1Algo/L1Field.cxx
@@ -16,7 +16,7 @@
// L1FieldValue methods
//
-bool L1FieldRegion::gkUseOriginalField = false;
+bool L1FieldRegion::fgkForceUseOfOriginalField = false;
//----------------------------------------------------------------------------------------------------------------------
//
@@ -173,24 +173,12 @@ void L1FieldRegion::CheckConsistency() const
// TODO: Any other checks? (S.Zharko)
}
-//----------------------------------------------------------------------------------------------------------------------
-// TODO: Should it be inline? (S.Zharko)
-L1FieldValue L1FieldRegion::Get(const fvec z)
-{
- fvec dz = z - z0;
- fvec dz2 = dz * dz;
- L1FieldValue B;
- B.x = cx0 + cx1 * dz + cx2 * dz2;
- B.y = cy0 + cy1 * dz + cy2 * dz2;
- B.z = cz0 + cz1 * dz + cz2 * dz2;
- return B;
-}
//----------------------------------------------------------------------------------------------------------------------
// TODO: Should it be inline? (S.Zharko)
-void L1FieldRegion::Get(const fvec x, const fvec y, const fvec z, fvec* B) const
+void L1FieldRegion::Get(const fvec& x, const fvec& y, const fvec& z, fvec* B) const
{
- if (gkUseOriginalField) {
+ if (fgkForceUseOfOriginalField || fUseOriginalField) {
for (size_t i = 0; i < fvec::size(); i++) {
double inPos[3] = {x[i], y[i], z[i]};
double outB[3];
@@ -201,26 +189,7 @@ void L1FieldRegion::Get(const fvec x, const fvec y, const fvec z, fvec* B) const
}
}
else {
- Get(z, B);
- }
-}
-
-//----------------------------------------------------------------------------------------------------------------------
-// TODO: Should it be inline? (S.Zharko)
-void L1FieldRegion::Get(const fvec z_, fvec* B) const
-{
- if (gkUseOriginalField) {
- for (size_t i = 0; i < fvec::size(); i++) {
- double inPos[3] = {0., 0., z_[i]};
- double outB[3];
- CbmKF::Instance()->GetMagneticField()->GetFieldValue(inPos, outB);
- B[0][i] = outB[0];
- B[1][i] = outB[1];
- B[2][i] = outB[2];
- }
- }
- else {
- fvec dz = z_ - z0;
+ fvec dz = z - z0;
fvec dz2 = dz * dz;
B[0] = cx0 + cx1 * dz + cx2 * dz2;
B[1] = cy0 + cy1 * dz + cy2 * dz2;
@@ -228,6 +197,7 @@ void L1FieldRegion::Get(const fvec z_, fvec* B) const
}
}
+
//----------------------------------------------------------------------------------------------------------------------
// TODO: Should it be inline? (S.Zharko)
void L1FieldRegion::Set(const L1FieldValue& b0, const fvec b0z, const L1FieldValue& b1, const fvec b1z,
diff --git a/reco/L1/L1Algo/L1Field.h b/reco/L1/L1Algo/L1Field.h
index 4d51fbc782..55efd242b9 100644
--- a/reco/L1/L1Algo/L1Field.h
+++ b/reco/L1/L1Algo/L1Field.h
@@ -104,21 +104,12 @@ public:
/// Consistency checker
void CheckConsistency() const;
- /// Gets the field vector at z
- // TODO: Probably we need a const specifier here, because the method does not change the fields
- L1FieldValue Get(const fvec z);
-
/// Gets the field vector and writes it into B pointer
/// \param x x-coordinate of the point to calculate the field
/// \param y y-coordinate of the point to calculate the field
/// \param z z-coordinate of the point to calculate the field
/// \param B pointer to the output fvec array of the magnetic field
- void Get(const fvec x, const fvec y, const fvec z, fvec* B) const;
-
- /// Gets the field vector and writes it into B pointer
- /// \param z_ z-coordinate of the point to calculate the field
- /// \param B pointer to the output fvec array of the magnetic field
- void Get(const fvec z_, fvec* B) const;
+ void Get(const fvec& x, const fvec& y, const fvec& z, fvec* B) const;
/// Interpolates the magnetic field by three nodal points and sets the result to this L1FieldRegion object
/// The result is a quadratic interpolation of the field as a function of z
@@ -170,8 +161,14 @@ public:
/// \param indentLevel number of indent characters in the output
std::string ToString(int indentLevel = 0) const;
+ /// Use original field instead of the field approximation
+ void SetUseOriginalField(bool v = true) { fUseOriginalField = v; }
+
+ /// Force using the original field
+ static void ForceUseOfOriginalField(bool v = true) { fgkForceUseOfOriginalField = v; }
+
public:
- static bool gkUseOriginalField;
+ static bool fgkForceUseOfOriginalField;
// TODO: Probably it's better to have arrays instead of separate fvec values? (S.Zharko)
// Bx(z) = cx0 + cx1*(z-z0) + cx2*(z-z0)^2
@@ -188,8 +185,11 @@ public:
fvec cz0 {0.f};
fvec cz1 {0.f};
fvec cz2 {0.f};
+
fvec z0 {0.f}; ///< z-coordinate of the field region central point
+ bool fUseOriginalField {false};
+
/// Serialization function
friend class boost::serialization::access;
template<class Archive>
diff --git a/reco/L1/L1Algo/L1InitManager.h b/reco/L1/L1Algo/L1InitManager.h
index 27820c2982..aff249d996 100644
--- a/reco/L1/L1Algo/L1InitManager.h
+++ b/reco/L1/L1Algo/L1InitManager.h
@@ -252,6 +252,9 @@ public:
/// Ignore hit search areas
void DevSetIgnoreHitSearchAreas(bool value = true) { fParameters.fDevIsIgnoreHitSearchAreas = value; }
+ /// Force use of the original field (not approximated)
+ void DevSetUseOfOriginalField(bool value = true) { fParameters.fDevIsUseOfOriginalField = value; }
+
/// Start singlets fit at the target
void DevSetFitSingletsFromTarget(bool value = true) { fParameters.fDevIsFitSingletsFromTarget = value; }
diff --git a/reco/L1/L1Algo/L1Parameters.h b/reco/L1/L1Algo/L1Parameters.h
index 7440474ae0..6111eac003 100644
--- a/reco/L1/L1Algo/L1Parameters.h
+++ b/reco/L1/L1Algo/L1Parameters.h
@@ -189,6 +189,9 @@ public:
/// Are the hit search areas ignored
bool DevIsIgnoreHitSearchAreas() const { return fDevIsIgnoreHitSearchAreas; }
+ /// Is original field must be used instead of the approximated one
+ bool DevIsUseOfOriginalField() const { return fDevIsUseOfOriginalField; }
+
/// Is singlets fit starts at the target
bool DevIsFitSingletsFromTarget() const { return fDevIsFitSingletsFromTarget; }
@@ -250,6 +253,7 @@ private:
// ***************************
bool fDevIsIgnoreHitSearchAreas {false}; ///< Process all hits on the station ignoring hit search area
+ bool fDevIsUseOfOriginalField {false}; ///< Force use of original field
bool fDevIsFitSingletsFromTarget {false}; ///< Fit singlet starting from the target with the KF
bool fDevIsMatchDoubletsViaMc {false}; ///< Flag to match doublets using MC information
bool fDevIsMatchTripletsViaMc {false}; ///< Flag to match triplets using Mc information
--
GitLab