-
Sergey Gorbunov authoredSergey Gorbunov authored
CaField.cxx 11.99 KiB
/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
SPDX-License-Identifier: GPL-3.0-only
Authors: Sergey Gorbunov, Sergei Zharko [committer] */
#include "CaField.h"
#include "KfUtils.h"
#include <iomanip>
#include <iostream>
#include <sstream>
using namespace cbm::algo::ca;
//
// FieldValue methods
//
template<typename DataT>
bool FieldRegion<DataT>::fgForceUseOfOriginalField = false;
template<typename DataT>
std::function<std::tuple<double, double, double>(double x, double y, double z)>
FieldRegion<DataT>::fgOriginalField([](double, double, double) {
assert(false && "ca::FieldRegion: The original field is not set, use cbm::ca::tools::SetCbmField()");
return std::tuple(kfdefs::Undef<double>, kfdefs::Undef<double>, kfdefs::Undef<double>);
});
template<typename DataT>
bool FieldRegion<DataT>::fgIsZeroOriginalField = false;
//----------------------------------------------------------------------------------------------------------------------
//
template<typename DataT>
void FieldValue<DataT>::CheckConsistency() const
{
/* Check SIMD data vectors for consistent initialization */
kfutils::CheckSimdVectorEquality(x, "FieldValue::x");
kfutils::CheckSimdVectorEquality(y, "FieldValue::y");
kfutils::CheckSimdVectorEquality(z, "FieldValue::z");
// TODO: Any other checks? (S.Zharko)
}
//----------------------------------------------------------------------------------------------------------------------
// TODO:
template<typename DataT>
std::string FieldValue<DataT>::ToString(int indentLevel) const
{
std::stringstream aStream{};
constexpr char indentChar = '\t';
std::string indent(indentLevel, indentChar);
aStream << indent << "Bx [kG]: " << std::setw(12) << std::setfill(' ') << kfutils::simd::Cast<DataT, float>(x)
<< '\n';
aStream << indent << "By [kG]: " << std::setw(12) << std::setfill(' ') << kfutils::simd::Cast<DataT, float>(y)
<< '\n';
aStream << indent << "Bz [kG]: " << std::setw(12) << std::setfill(' ') << kfutils::simd::Cast<DataT, float>(z);
return aStream.str();
}
namespace cbm::algo::ca
{
template class FieldValue<fvec>;
template class FieldValue<float>;
template class FieldValue<double>;
} // namespace cbm::algo::ca
//
// 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;
B.x = 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;
B.y = 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;
B.z = 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
//
//----------------------------------------------------------------------------------------------------------------------
//
template<typename DataT>
FieldRegion<DataT>::FieldRegion(float reg[10]) noexcept
: cx0(reg[0])
, cx1(reg[1])
, cx2(reg[2])
, cy0(reg[3])
, cy1(reg[4])
, cy2(reg[5])
, cz0(reg[6])
, cz1(reg[7])
, cz2(reg[8])
, z0(reg[9])
{
}
//----------------------------------------------------------------------------------------------------------------------
//
template<typename DataT>
void FieldRegion<DataT>::CheckConsistency() const
{
/* Check SIMD data vectors for consistent initialization */
kfutils::CheckSimdVectorEquality(cx0, "FieldRegion::cx0");
kfutils::CheckSimdVectorEquality(cx1, "FieldRegion::cx1");
kfutils::CheckSimdVectorEquality(cx2, "FieldRegion::cx2");
kfutils::CheckSimdVectorEquality(cy0, "FieldRegion::cy0");
kfutils::CheckSimdVectorEquality(cy1, "FieldRegion::cy1");
kfutils::CheckSimdVectorEquality(cy2, "FieldRegion::cy2");
kfutils::CheckSimdVectorEquality(cz0, "FieldRegion::cz0");
kfutils::CheckSimdVectorEquality(cz1, "FieldRegion::cz1");
kfutils::CheckSimdVectorEquality(cz2, "FieldRegion::cz2"); kfutils::CheckSimdVectorEquality(z0, "FieldRegion::z0");
// TODO: Any other checks? (S.Zharko)
}
//----------------------------------------------------------------------------------------------------------------------
// TODO: Should it be inline? (S.Zharko)
template<typename DataT>
FieldValue<DataT> FieldRegion<DataT>::Get(const DataT& x, const DataT& y, const DataT& z) const
{
FieldValue<DataT> B;
if (fgForceUseOfOriginalField || fUseOriginalField) {
// TODO: ? Why the simd words are set indvidually?
for (size_t i = 0; i < kfutils::simd::Size<DataT>(); i++) {
auto [bx, by, bz] =
fgOriginalField(kfutils::simd::Cast<DataT, double>(x, i), kfutils::simd::Cast<DataT, double>(y, i),
kfutils::simd::Cast<DataT, double>(z, i));
B.SetSimdEntry(bx, by, bz, i);
}
}
else {
DataT dz = z - z0;
DataT dz2 = dz * dz;
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)
template<typename DataT>
void FieldRegion<DataT>::Set(const FieldValue<DataT>& b0, const DataT b0z, const FieldValue<DataT>& b1, const DataT b1z,
const FieldValue<DataT>& b2, const DataT b2z)
{
fIsZeroField = (b0.IsZeroField() && b1.IsZeroField() && b2.IsZeroField());
z0 = b0z;
DataT dz1 = b1z - b0z;
DataT dz2 = b2z - b0z;
DataT det = kfutils::simd::One<DataT>() / (dz1 * dz2 * (b2z - b1z));
DataT w21 = -dz2 * det;
DataT w22 = dz1 * det;
DataT w11 = -dz2 * w21;
DataT w12 = -dz1 * w22;
DataT db1 = b1.x - b0.x;
DataT db2 = b2.x - b0.x;
cx0 = b0.x;
cx1 = db1 * w11 + db2 * w12;
cx2 = db1 * w21 + db2 * w22;
db1 = b1.y - b0.y;
db2 = b2.y - b0.y;
cy0 = b0.y;
cy1 = db1 * w11 + db2 * w12;
cy2 = db1 * w21 + db2 * w22;
db1 = b1.z - b0.z;
db2 = b2.z - b0.z;
cz0 = b0.z;
cz1 = db1 * w11 + db2 * w12;
cz2 = db1 * w21 + db2 * w22;
}
//----------------------------------------------------------------------------------------------------------------------
// TODO: Should it be inline? (S.Zharko)//template<typename DataT>
//void FieldRegion<DataT>::Set(const FieldValue<DataT>& b0, const DataT b0z, const FieldValue<DataT>& b1, const DataT b1z)
//{
// fIsZeroField = (b0.IsZeroField() && b1.IsZeroField());
//
// z0 = b0z;
// DataT dzi = kfutils::simd::One<DataT>() / (b1z - b0z);
// cx0 = b0.x;
// cy0 = b0.y;
// cz0 = b0.z;
// cx1 = (b1.x - b0.x) * dzi;
// cy1 = (b1.y - b0.y) * dzi;
// cz1 = (b1.z - b0.z) * dzi;
// cx2 = 0.f;
// cy2 = 0.f;
// cz2 = 0.f;
//}
//----------------------------------------------------------------------------------------------------------------------
// TODO: Should it be inline? (S.Zharko)
template<typename DataT>
void FieldRegion<DataT>::Shift(DataT z)
{
DataT dz = z - z0;
DataT cx2dz = cx2 * dz;
DataT cy2dz = cy2 * dz;
DataT cz2dz = cz2 * dz;
z0 = z;
cx0 += (cx1 + cx2dz) * dz;
cy0 += (cy1 + cy2dz) * dz;
cz0 += (cz1 + cz2dz) * dz;
cx1 += cx2dz + cx2dz;
cy1 += cy2dz + cy2dz;
cz1 += cz2dz + cz2dz;
}
//----------------------------------------------------------------------------------------------------------------------
//
template<typename DataT>
std::string FieldRegion<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 << "Bx(z) [kG] components: " << '\n';
aStream << indent << indentChar << "cx0: " << cx0 << '\n';
aStream << indent << indentChar << "cx1: " << cx1 << '\n';
aStream << indent << indentChar << "cx2: " << cx2 << '\n';
aStream << indent << "By(z) [kG] components: " << '\n';
aStream << indent << indentChar << "cy0: " << cy0 << '\n';
aStream << indent << indentChar << "cy1: " << cy1 << '\n';
aStream << indent << indentChar << "cy2: " << cy2 << '\n';
aStream << indent << "Bz(z) [kG] components: " << '\n';
aStream << indent << indentChar << "cz0: " << cz0 << '\n';
aStream << indent << indentChar << "cz1: " << cz1 << '\n';
aStream << indent << indentChar << "cz2: " << cz2 << '\n';
aStream << indent << "z0 [cm]: " << z0;
return aStream.str();
}
namespace cbm::algo::ca
{
template class FieldRegion<fvec>;
template class FieldRegion<float>;
template class FieldRegion<double>;
} // namespace cbm::algo::ca