diff --git a/algo/ca/core/pars/CaField.cxx b/algo/ca/core/pars/CaField.cxx
index b2666bd454bd0fbdb2692cef8949267553037697..676a53135bedacfde382cdc1c72422105e0aef49 100644
--- a/algo/ca/core/pars/CaField.cxx
+++ b/algo/ca/core/pars/CaField.cxx
@@ -12,9 +12,7 @@
using namespace cbm::algo::ca;
-//
-// FieldValue methods
-//
+
template<typename DataT>
bool FieldRegion<DataT>::fgForceUseOfOriginalField = false;
@@ -29,117 +27,6 @@ template<typename DataT>
bool FieldRegion<DataT>::fgIsZeroOriginalField = false;
-//
-// FieldSlice methods
-//
-
-//----------------------------------------------------------------------------------------------------------------------
-//
-template<typename DataT>
-FieldSlice<DataT>::FieldSlice()
-{
- for (int i = 0; i < constants::size::MaxNFieldApproxCoefficients; ++i) {
- cx[i] = constants::Undef<fscal>;
- cy[i] = constants::Undef<fscal>;
- cz[i] = constants::Undef<fscal>;
- }
- z = constants::Undef<fscal>;
-}
-
-//----------------------------------------------------------------------------------------------------------------------
-//
-template<typename DataT>
-void FieldSlice<DataT>::CheckConsistency() const
-{
- /* Check SIMD data vectors for consistent initialization */
- for (int i = 0; i < constants::size::MaxNFieldApproxCoefficients; ++i) {
- kfutils::CheckSimdVectorEquality(cx[i], "FieldSlice: cx");
- kfutils::CheckSimdVectorEquality(cy[i], "FieldSlice: cy");
- kfutils::CheckSimdVectorEquality(cz[i], "FieldSlice: cz");
- }
- kfutils::CheckSimdVectorEquality(z, "FieldSlice: z");
-}
-
-
-//----------------------------------------------------------------------------------------------------------------------
-// TODO: Should it be inline? (S.Zharko)
-template<typename DataT>
-FieldValue<DataT> FieldSlice<DataT>::GetFieldValue(const DataT& x, const DataT& y) const
-{
- DataT x2 = x * x;
- DataT y2 = y * y;
- DataT xy = x * y;
-
- DataT x3 = x2 * x;
- DataT y3 = y2 * y;
- DataT xy2 = x * y2;
- DataT x2y = x2 * y;
-
- DataT x4 = x3 * x;
- DataT y4 = y3 * y;
- DataT xy3 = x * y3;
- DataT x2y2 = x2 * y2;
- DataT x3y = x3 * y;
-
- DataT x5 = x4 * x;
- DataT y5 = y4 * y;
- DataT xy4 = x * y4;
- DataT x2y3 = x2 * y3;
- DataT x3y2 = x3 * y2;
- DataT x4y = x4 * y;
-
- FieldValue<DataT> B(
- // Bx
- cx[0] + cx[1] * x + cx[2] * y + cx[3] * x2 + cx[4] * xy + cx[5] * y2 + cx[6] * x3 + cx[7] * x2y + cx[8] * xy2
- + cx[9] * y3 + cx[10] * x4 + cx[11] * x3y + cx[12] * x2y2 + cx[13] * xy3 + cx[14] * y4 + cx[15] * x5
- + cx[16] * x4y + cx[17] * x3y2 + cx[18] * x2y3 + cx[19] * xy4 + cx[20] * y5,
- // By
- cy[0] + cy[1] * x + cy[2] * y + cy[3] * x2 + cy[4] * xy + cy[5] * y2 + cy[6] * x3 + cy[7] * x2y + cy[8] * xy2
- + cy[9] * y3 + cy[10] * x4 + cy[11] * x3y + cy[12] * x2y2 + cy[13] * xy3 + cy[14] * y4 + cy[15] * x5
- + cy[16] * x4y + cy[17] * x3y2 + cy[18] * x2y3 + cy[19] * xy4 + cy[20] * y5,
- // Bz
- cz[0] + cz[1] * x + cz[2] * y + cz[3] * x2 + cz[4] * xy + cz[5] * y2 + cz[6] * x3 + cz[7] * x2y + cz[8] * xy2
- + cz[9] * y3 + cz[10] * x4 + cz[11] * x3y + cz[12] * x2y2 + cz[13] * xy3 + cz[14] * y4 + cz[15] * x5
- + cz[16] * x4y + cz[17] * x3y2 + cz[18] * x2y3 + cz[19] * xy4 + cz[20] * y5);
-
- return B;
-}
-
-template<typename DataT>
-FieldValue<DataT> FieldSlice<DataT>::GetFieldValueForLine(const TrackParamBase<DataT>& t) const
-{
- DataT dz = z - t.GetZ();
- return GetFieldValue(t.GetX() + t.GetTx() * dz, t.GetY() + t.GetTy() * dz);
-}
-
-
-//----------------------------------------------------------------------------------------------------------------------
-//
-template<typename DataT>
-std::string FieldSlice<DataT>::ToString(int indentLevel) const
-{
- std::stringstream aStream{};
- // TODO: possibly it is better to place the indentChar into L1Parameters (S.Zharko)
- constexpr char indentChar = '\t';
- std::string indent(indentLevel, indentChar);
- aStream << indent << "idx CX CY CZ";
- for (int i = 0; i < constants::size::MaxNFieldApproxCoefficients; ++i) {
- aStream << '\n' << indent;
- aStream << std::setw(3) << std::setfill(' ') << i << ' ';
- aStream << std::setw(12) << std::setfill(' ') << kfutils::simd::Cast<DataT, float>(cx[i], 0) << ' ';
- aStream << std::setw(12) << std::setfill(' ') << kfutils::simd::Cast<DataT, float>(cy[i], 0) << ' ';
- aStream << std::setw(12) << std::setfill(' ') << kfutils::simd::Cast<DataT, float>(cz[i], 0);
- }
- return aStream.str();
-}
-
-namespace cbm::algo::ca
-{
- template class FieldSlice<fvec>;
- template class FieldSlice<float>;
- template class FieldSlice<double>;
-} // namespace cbm::algo::ca
-
//
// FieldRegion methdos
//
diff --git a/algo/ca/core/pars/CaField.h b/algo/ca/core/pars/CaField.h
index a12ae18e5b5010ce506ff44e6cd993ab5d8aa055..023eed73858fec5ceadba4c81e17b6738865ed36 100644
--- a/algo/ca/core/pars/CaField.h
+++ b/algo/ca/core/pars/CaField.h
@@ -7,6 +7,7 @@
#include "CaDefs.h"
#include "CaSimd.h"
+#include "KfFieldValue.h"
#include "KfTrackParam.h"
#include <boost/serialization/access.hpp>
@@ -16,68 +17,9 @@
namespace cbm::algo::ca
{
+ using cbm::algo::kf::FieldSlice;
using cbm::algo::kf::FieldValue;
- /// Class represents a set of magnetic field approximation coefficients
- ///
- // TODO: Crosscheck the default content (S.Zharko)
- template<typename DataT>
- class FieldSlice {
- public:
- /// Default constructor
- FieldSlice();
-
- /// \brief Copy constructor with type conversion
- template<typename DataIn>
- FieldSlice(const FieldSlice<DataIn>& other) : z(kfutils::simd::Cast<DataIn, DataT>(other.z))
- {
- for (size_t i = 0; i < kMaxNFieldApproxCoefficients; i++) {
- cx[i] = kfutils::simd::Cast<DataIn, DataT>(other.cx[i]);
- cy[i] = kfutils::simd::Cast<DataIn, DataT>(other.cy[i]);
- cz[i] = kfutils::simd::Cast<DataIn, DataT>(other.cz[i]);
- }
- }
-
- /// Consistency checker
- void CheckConsistency() const;
-
- /// Gets field value from (x, y) fvec point
- /// \param x x-coordinate of input
- /// \param y y-coordinate of input
- /// \return B the FieldValue output
- FieldValue<DataT> GetFieldValue(const DataT& x, const DataT& y) const;
-
- /// Gets field value for the intersection with a straight track
- /// \param t straight track
- /// \return B the FieldValue output
- FieldValue<DataT> GetFieldValueForLine(const TrackParamBase<DataT>& t) const;
-
- /// String representation of class contents
- /// \param indentLevel number of indent characters in the output
- std::string ToString(int indentLevel = 0) const;
-
- public:
- // NOTE: We don't use an initialization of arrays here because we cannot be sure
- // if the underlying type (fvec) has a default constructor, but
- // we are sure, that it can be initialized with a float. (S.Zharko)
- static constexpr auto kMaxNFieldApproxCoefficients = constants::size::MaxNFieldApproxCoefficients;
- DataT cx[kMaxNFieldApproxCoefficients]; ///< Polynomial coefficients for x-component of the field value
- DataT cy[kMaxNFieldApproxCoefficients]; ///< Polynomial coefficients for y-component of the field value
- DataT cz[kMaxNFieldApproxCoefficients]; ///< Polynomial coefficients for z-component of the field value
-
- DataT z{constants::Undef<float>}; ///< z coordinate of the slice
-
- /// Serialization function
- friend class boost::serialization::access;
- template<class Archive>
- void serialize(Archive& ar, const unsigned int)
- {
- ar& cx;
- ar& cy;
- ar& cz;
- ar& z;
- }
- } _fvecalignment;
template<typename DataT>
class FieldRegion {
diff --git a/algo/ca/core/pars/CaStationInitializer.cxx b/algo/ca/core/pars/CaStationInitializer.cxx
index 54ae533c1c0e0fc9e8ae3be87be2fc0fdc1e8056..3052dc9662b47f558ccdb752e962e8d11a13f09e 100644
--- a/algo/ca/core/pars/CaStationInitializer.cxx
+++ b/algo/ca/core/pars/CaStationInitializer.cxx
@@ -108,76 +108,17 @@ void StationInitializer::SetFieldFunction(
if (!fInitController.GetFlag(EInitKey::kYmax)) {
LOG(fatal) << "Attempt to set magnetic field slice before Ymax size of the station";
}
- // TODO: Change names of variables according to convention (S.Zh.)
- constexpr int M = constants::size::MaxFieldApproxPolynomialOrder;
- constexpr int N = constants::size::MaxNFieldApproxCoefficients;
- constexpr int D = 3; ///> number of dimensions
-
- // SLE initialization
- double A[N][N + D] = {}; // augmented matrix
- double dx = (fXmax / N / 2 < 1.) ? fXmax / N / 4. : 1.;
- double dy = (fYmax / N / 2 < 1.) ? fYmax / N / 4. : 1.;
-
- for (double x = -fXmax; x <= fXmax; x += dx) {
- for (double y = -fYmax; y <= fYmax; y += dy) {
- double r = sqrt(fabs(x * x / fXmax / fXmax + y / fYmax * y / fYmax));
- if (r > 1.) {
- continue;
- }
- double p[D] = {x, y, fZref};
- double B[D] = {};
- getFieldValue(p, B);
-
- double m[N] = {1};
- m[0] = 1;
- for (int i = 1; i <= M; ++i) {
- int k = (i - 1) * i / 2;
- int l = i * (i + 1) / 2;
- for (int j = 0; j < i; ++j) {
- m[l + j] = x * m[k + j];
- }
- m[l + i] = y * m[k + i - 1];
- }
-
- double w = 1. / (r * r + 1);
- for (int i = 0; i < N; ++i) {
- // fill the left part of the matrix
- for (int j = 0; j < N; ++j) {
- A[i][j] += w * m[i] * m[j];
- }
- // fill the right part of the matrix
- for (int j = 0; j < D; ++j) {
- A[i][N + j] += w * B[j] * m[i];
- }
- }
- }
- }
- // SLE solution (Gaussian elimination)
- //
- for (int kCol = 0; kCol < N - 1; ++kCol) {
- for (int jRow = kCol + 1; jRow < N; ++jRow) {
- double factor = A[jRow][kCol] / A[kCol][kCol];
- for (int iCol = kCol; iCol < N + D; ++iCol) {
- A[jRow][iCol] -= factor * A[kCol][iCol];
- }
- }
- }
- for (int kCol = N - 1; kCol > 0; --kCol) {
- for (int jRow = kCol - 1; jRow >= 0; --jRow) {
- double factor = A[jRow][kCol] / A[kCol][kCol];
- for (int iCol = kCol; iCol < N + D; ++iCol) {
- A[jRow][iCol] -= factor * A[kCol][iCol];
- }
- }
- }
+ /// \brief Magnetic field function type
+ /// Signature: tuple<Bx, By, Bz>(x, y, z);
+ auto glambda = [&](double x, double y, double z) -> std::tuple<double, double, double> {
+ double xyz[3] = {x, y, z};
+ double B[3] = {};
+ getFieldValue(xyz, B);
+ return std::tuple<double, double, double>(B[0], B[1], B[2]);
+ };
- for (int j = 0; j < N; ++j) {
- fStation.fieldSlice.cx[j] = A[j][N + 0] / A[j][j];
- fStation.fieldSlice.cy[j] = A[j][N + 1] / A[j][j];
- fStation.fieldSlice.cz[j] = A[j][N + 2] / A[j][j];
- }
- fStation.fieldSlice.z = fZref;
+ fStation.fieldSlice = FieldSlice<fvec>(glambda, fXmax, fYmax, fZref);
fInitController.SetFlag(EInitKey::kFieldSlice);
}
diff --git a/algo/ca/core/pars/CaStationInitializer.h b/algo/ca/core/pars/CaStationInitializer.h
index 5cf515bc8ff7f51ac5fbb060cc8434ff78c03dc1..6f328f2f7eae2192f488c1230cb7c194c2558824 100644
--- a/algo/ca/core/pars/CaStationInitializer.h
+++ b/algo/ca/core/pars/CaStationInitializer.h
@@ -75,25 +75,6 @@ namespace cbm::algo::ca
/// \brief Gets detector ID
EDetectorID GetDetectorID() const { return fDetectorID; }
- /// \brief Gets a coefficient with idx for the field x-component approximation
- fvec GetFieldSliceCx(int idx) const { return fStation.fieldSlice.cx[idx]; }
-
- /// \brief Gets a coefficient with idx for the field y-component approximation
- fvec GetFieldSliceCy(int idx) const { return fStation.fieldSlice.cy[idx]; }
-
- /// \brief Gets a coefficient with idx for the field z-component approximation
- fvec GetFieldSliceCz(int idx) const { return fStation.fieldSlice.cz[idx]; }
-
- /// \brief Gets array of the coefficients for the field x-component approximation
- const fvec* GetFieldSliceCx() const { return fStation.fieldSlice.cx; }
-
- /// \brief Gets array of the coefficients for the field y-component approximation
- const fvec* GetFieldSliceCy() const { return fStation.fieldSlice.cy; }
-
- /// \brief Gets array of the coefficients for the field z-component approximation
- const fvec* GetFieldSliceCz() const { return fStation.fieldSlice.cz; }
-
-
/// \brief Gets field status: 0 - station is outside the field, 1 - station is inside the field
int GetFieldStatus() const { return fStation.fieldStatus; }
diff --git a/algo/ca/core/tracking/CaTrackFitter.cxx b/algo/ca/core/tracking/CaTrackFitter.cxx
index 6b3b8407459195a32185744940b6aaf9d6eb1083..a2f3214b9b1a92a3f15b238c5f4627491c693e47 100644
--- a/algo/ca/core/tracking/CaTrackFitter.cxx
+++ b/algo/ca/core/tracking/CaTrackFitter.cxx
@@ -203,7 +203,7 @@ void TrackFitter::FitCaTracks(const ca::InputData& input, WindowData& wData)
for (int ih = nHitsTrack - 1; ih >= 0; ih--) {
const int ista = iSta[ih];
const ca::Station<fvec>& st = fParameters.GetStation(ista);
- By[ista][iVec] = st.fieldSlice.cy[0][0];
+ By[ista][iVec] = st.fieldSlice.GetFieldValue(0., 0.).GetBy()[0];
}
}
diff --git a/algo/kf/core/geo/KfFieldSlice.cxx b/algo/kf/core/geo/KfFieldSlice.cxx
index 8e8989ea1c8d0c212853abe841efafc92e7ae243..38d8bd9828301e01ae0859fb95816689faacb900 100644
--- a/algo/kf/core/geo/KfFieldSlice.cxx
+++ b/algo/kf/core/geo/KfFieldSlice.cxx
@@ -28,8 +28,15 @@ FieldSlice<T>::FieldSlice(const FieldFn_t& fieldFn, double xMax, double yMax, do
//
std::array<std::array<double, kNofCoeff + 3>, kNofCoeff> A = {{}};
- const double dx = (xMax < 2. * kNofCoeff) ? (xMax / kNofCoeff / 4.) : 1.;
- const double dy = (yMax < 2. * kNofCoeff) ? (yMax / kNofCoeff / 4.) : 1.;
+ double dx = xMax / kNofCoeff / 10.;
+ double dy = yMax / kNofCoeff / 10.;
+
+ if (dx > 1.) {
+ dx = 1.;
+ }
+ if (dy > 1.) {
+ dy = 1.;
+ }
// Point and field arrays to access the field
std::array<double, 3> field;
@@ -40,11 +47,6 @@ FieldSlice<T>::FieldSlice(const FieldFn_t& fieldFn, double xMax, double yMax, do
// Fill the augmented matrix
for (double x = -xMax; x <= xMax; x += dx) {
for (double y = -yMax; y <= yMax; y += dy) {
- const double r2 = x * x / (xMax * xMax) + y * y / (yMax * yMax);
- if (r2 > 1.) {
- continue;
- }
-
std::tie(field[0], field[1], field[2]) = fieldFn(x, y, zRef);
// Fill the monomial array
{
@@ -63,7 +65,7 @@ FieldSlice<T>::FieldSlice(const FieldFn_t& fieldFn, double xMax, double yMax, do
}
{
- double w = 1. / (r2 + 1.);
+ double w = 1.; // / (r2 + 1.);
for (int i = 0; i < kNofCoeff; ++i) {
// Fill the left part
for (int j = 0; j < kNofCoeff; ++j) {
@@ -143,6 +145,18 @@ std::string FieldSlice<T>::ToString(int indentLevel, int verbose) const
return msg.str();
}
+template<typename T>
+void FieldSlice<T>::CheckConsistency() const
+{
+ // Check SIMD data vectors for consistent initialization
+ for (int i = 0; i < kNofCoeff; ++i) {
+ utils::CheckSimdVectorEquality(fBx[i], "FieldSlice: cx");
+ utils::CheckSimdVectorEquality(fBy[i], "FieldSlice: cy");
+ utils::CheckSimdVectorEquality(fBz[i], "FieldSlice: cz");
+ }
+ utils::CheckSimdVectorEquality(fZref, "FieldSlice: z");
+}
+
namespace cbm::algo::kf
{
template class FieldSlice<float>;
diff --git a/algo/kf/core/geo/KfFieldSlice.h b/algo/kf/core/geo/KfFieldSlice.h
index 6be46bad7ce19241551da0f348dea12fecc890e7..34a6e201abdeefb3c5e0b9dd75bebeac628ffed2 100644
--- a/algo/kf/core/geo/KfFieldSlice.h
+++ b/algo/kf/core/geo/KfFieldSlice.h
@@ -73,11 +73,9 @@ namespace cbm::algo::kf
constexpr FieldValue<T> GetFieldValue(const T& x, const T& y) const
{
using math::Horner;
- /* clang-format off */
- return FieldValue<T>(Horner<kPolDegree>(fBx.cbegin(), x, y),
- Horner<kPolDegree>(fBy.cbegin(), x, y),
+ return FieldValue<T>(Horner<kPolDegree>(fBx.cbegin(), x, y), //
+ Horner<kPolDegree>(fBy.cbegin(), x, y), //
Horner<kPolDegree>(fBz.cbegin(), x, y));
- /* clang-format on */
}
/// \brief Gets field value for the intersection with a straight track
@@ -97,6 +95,9 @@ namespace cbm::algo::kf
/// \param verbose Verbosity level
std::string ToString(int indentLevel = 0, int verbose = 1) const;
+ /// Consistency checker
+ void CheckConsistency() const;
+
private:
/// \brief Serialization function
friend class boost::serialization::access;
diff --git a/reco/L1/CbmL1Performance.cxx b/reco/L1/CbmL1Performance.cxx
index 6a43c542099cb4d0b6653bdea1c0efb17b2b19d4..92ac7e05a2b809388037707af92b4a064cd3ec10 100644
--- a/reco/L1/CbmL1Performance.cxx
+++ b/reco/L1/CbmL1Performance.cxx
@@ -1572,7 +1572,7 @@ void CbmL1::FieldApproxCheck()
{
TDirectory* curr = gDirectory;
TFile* currentFile = gFile;
- TFile* fout = new TFile("FieldApprox.root", "RECREATE");
+ TFile* fout = new TFile("CaFieldApproxQa.root", "RECREATE");
fout->cd();
assert(FairRunAna::Instance());
@@ -1590,19 +1590,56 @@ void CbmL1::FieldApproxCheck()
int NbinsX = 101;
int NbinsY = 101;
- TH2F* stB = new TH2F(Form("station_%i_dB", ist), Form("station %i, dB, z = %0.f cm", ist, z), //
- NbinsX, -Xmax, Xmax, NbinsY, -Ymax, Ymax);
- TH2F* stBx = new TH2F(Form("station_%i_dBx", ist), Form("station %i, dBx, z = %0.f cm", ist, z), //
- NbinsX, -Xmax, Xmax, NbinsY, -Ymax, Ymax);
- TH2F* stBy = new TH2F(Form("station_%i_dBy", ist), Form("station %i, dBy, z = %0.f cm", ist, z), //
- NbinsX, -Xmax, Xmax, NbinsY, -Ymax, Ymax);
- TH2F* stBz = new TH2F(Form("station_%i_dBz", ist), Form("station %i, dBz, z = %0.f cm", ist, z), //
- NbinsX, -Xmax, Xmax, NbinsY, -Ymax, Ymax);
-
- for (int i = 0; i < stB->GetXaxis()->GetNbins(); i++) {
- double x = stB->GetXaxis()->GetBinCenter(i);
- for (int j = 0; j < stB->GetYaxis()->GetNbins(); j++) {
- double y = stB->GetYaxis()->GetBinCenter(j);
+ TH2F* hdB[4] = {
+ new TH2F(Form("station_%i_dB", ist), Form("station %i, dB, z = %0.f cm", ist, z), //
+ NbinsX, -Xmax, Xmax, NbinsY, -Ymax, Ymax),
+ new TH2F(Form("station_%i_dBx", ist), Form("station %i, dBx, z = %0.f cm", ist, z), //
+ NbinsX, -Xmax, Xmax, NbinsY, -Ymax, Ymax),
+ new TH2F(Form("station_%i_dBy", ist), Form("station %i, dBy, z = %0.f cm", ist, z), //
+ NbinsX, -Xmax, Xmax, NbinsY, -Ymax, Ymax),
+ new TH2F(Form("station_%i_dBz", ist), Form("station %i, dBz, z = %0.f cm", ist, z), //
+ NbinsX, -Xmax, Xmax, NbinsY, -Ymax, Ymax) //
+ };
+
+ TH2F* hB[4] = {
+ new TH2F(Form("station_%i_B", ist), Form("station %i, B, z = %0.f cm", ist, z), //
+ NbinsX, -Xmax, Xmax, NbinsY, -Ymax, Ymax),
+ new TH2F(Form("station_%i_Bx", ist), Form("station %i, Bx, z = %0.f cm", ist, z), //
+ NbinsX, -Xmax, Xmax, NbinsY, -Ymax, Ymax),
+ new TH2F(Form("station_%i_By", ist), Form("station %i, By, z = %0.f cm", ist, z), //
+ NbinsX, -Xmax, Xmax, NbinsY, -Ymax, Ymax),
+ new TH2F(Form("station_%i_Bz", ist), Form("station %i, Bz, z = %0.f cm", ist, z), //
+ NbinsX, -Xmax, Xmax, NbinsY, -Ymax, Ymax) //
+ };
+
+ for (int i = 0; i < 4; i++) {
+ hdB[i]->GetXaxis()->SetTitle("X, cm");
+ hdB[i]->GetYaxis()->SetTitle("Y, cm");
+ hdB[i]->GetXaxis()->SetTitleOffset(1);
+ hdB[i]->GetYaxis()->SetTitleOffset(1);
+ hdB[i]->GetZaxis()->SetTitleOffset(1.3);
+ hB[i]->GetXaxis()->SetTitle("X, cm");
+ hB[i]->GetYaxis()->SetTitle("Y, cm");
+ hB[i]->GetXaxis()->SetTitleOffset(1);
+ hB[i]->GetYaxis()->SetTitleOffset(1);
+ hB[i]->GetZaxis()->SetTitleOffset(1.3);
+ }
+
+ hdB[0]->GetZaxis()->SetTitle("B_map - B_L1, kGauss");
+ hdB[1]->GetZaxis()->SetTitle("Bx_map - Bx_L1, kGauss");
+ hdB[2]->GetZaxis()->SetTitle("By_map - By_L1, kGauss");
+ hdB[3]->GetZaxis()->SetTitle("Bz_map - Bz_L1, kGauss");
+
+ hdB[0]->GetZaxis()->SetTitle("B_map, kGauss");
+ hdB[1]->GetZaxis()->SetTitle("Bx_map, kGauss");
+ hdB[2]->GetZaxis()->SetTitle("By_map, kGauss");
+ hdB[3]->GetZaxis()->SetTitle("Bz_map, kGauss");
+
+
+ for (int i = 1; i <= hdB[0]->GetXaxis()->GetNbins(); i++) {
+ double x = hdB[0]->GetXaxis()->GetBinCenter(i);
+ for (int j = 1; j <= hdB[0]->GetYaxis()->GetNbins(); j++) {
+ double y = hdB[0]->GetYaxis()->GetBinCenter(j);
double r[] = {x, y, z};
double B[] = {0., 0., 0.};
MF->GetFieldValue(r, B);
@@ -1611,45 +1648,22 @@ void CbmL1::FieldApproxCheck()
ca::FieldValue<fvec> B_L1 = st.fieldSlice.GetFieldValue(x, y);
- stB->SetBinContent(i, j, (bbb - B_L1.GetAbs()[0]));
- stBx->SetBinContent(i, j, (B[0] - B_L1.GetBx()[0]));
- stBy->SetBinContent(i, j, (B[1] - B_L1.GetBy()[0]));
- stBz->SetBinContent(i, j, (B[2] - B_L1.GetBz()[0]));
+ hdB[0]->SetBinContent(i, j, (bbb - B_L1.GetAbs()[0]));
+ hdB[1]->SetBinContent(i, j, (B[0] - B_L1.GetBx()[0]));
+ hdB[2]->SetBinContent(i, j, (B[1] - B_L1.GetBy()[0]));
+ hdB[3]->SetBinContent(i, j, (B[2] - B_L1.GetBz()[0]));
+
+ hB[0]->SetBinContent(i, j, bbb);
+ hB[1]->SetBinContent(i, j, B[0]);
+ hB[2]->SetBinContent(i, j, B[1]);
+ hB[3]->SetBinContent(i, j, B[2]);
}
}
- stB->GetXaxis()->SetTitle("X, cm");
- stB->GetYaxis()->SetTitle("Y, cm");
- stB->GetXaxis()->SetTitleOffset(1);
- stB->GetYaxis()->SetTitleOffset(1);
- stB->GetZaxis()->SetTitle("B_map - B_L1, kGauss");
- stB->GetZaxis()->SetTitleOffset(1.3);
-
- stBx->GetXaxis()->SetTitle("X, cm");
- stBx->GetYaxis()->SetTitle("Y, cm");
- stBx->GetXaxis()->SetTitleOffset(1);
- stBx->GetYaxis()->SetTitleOffset(1);
- stBx->GetZaxis()->SetTitle("Bx_map - Bx_L1, kGauss");
- stBx->GetZaxis()->SetTitleOffset(1.3);
-
- stBy->GetXaxis()->SetTitle("X, cm");
- stBy->GetYaxis()->SetTitle("Y, cm");
- stBy->GetXaxis()->SetTitleOffset(1);
- stBy->GetYaxis()->SetTitleOffset(1);
- stBy->GetZaxis()->SetTitle("By_map - By_L1, kGauss");
- stBy->GetZaxis()->SetTitleOffset(1.3);
-
- stBz->GetXaxis()->SetTitle("X, cm");
- stBz->GetYaxis()->SetTitle("Y, cm");
- stBz->GetXaxis()->SetTitleOffset(1);
- stBz->GetYaxis()->SetTitleOffset(1);
- stBz->GetZaxis()->SetTitle("Bz_map - Bz_L1, kGauss");
- stBz->GetZaxis()->SetTitleOffset(1.3);
-
- stB->Write();
- stBx->Write();
- stBy->Write();
- stBz->Write();
+ for (int i = 0; i < 4; i++) {
+ hdB[i]->Write();
+ hB[i]->Write();
+ }
} // for ista