diff --git a/MagneticField/ParametrizedEngine/interface/ParabolicParametrizedMagneticField.h b/MagneticField/ParametrizedEngine/interface/ParabolicParametrizedMagneticField.h deleted file mode 100644 index 1107e6ee4cd85..0000000000000 --- a/MagneticField/ParametrizedEngine/interface/ParabolicParametrizedMagneticField.h +++ /dev/null @@ -1,47 +0,0 @@ -/** - Description: Utility function to calculate the Magnetic Field on the GPU. The Vec3 argument of the functions must support access to its components via (), note that e.g. Eigen::Matrix provides such an interface. -*/ - -#ifndef MagneticField_ParametrizedEngine_interface_ParabolicParametrizedMagneticField_h -#define MagneticField_ParametrizedEngine_interface_ParabolicParametrizedMagneticField_h - -namespace magneticFieldParabolicPortable { - - struct Parameters { - // These parameters are the best fit of 3.8T to the OAEParametrizedMagneticField parametrization. - // See MagneticField/ParametrizedEngine/src/ParabolicParametrizedMagneticField.cc - static constexpr float c1 = 3.8114; - static constexpr float b0 = -3.94991e-06; - static constexpr float b1 = 7.53701e-06; - static constexpr float a = 2.43878e-11; - static constexpr float max_radius2 = 13225.f; // tracker radius - static constexpr float max_z = 280.f; // tracker z - }; - - template - constexpr float Kr(Vec3 const& vec) { - return Parameters::a * (vec(0) * vec(0) + vec(1) * vec(1)) + 1.f; - } - - template - constexpr float B0Z(Vec3 const& vec) { - return Parameters::b0 * vec(2) * vec(2) + Parameters::b1 * vec(2) + Parameters::c1; - } - - template - constexpr bool isValid(Vec3 const& vec) { - return ((vec(0) * vec(0) + vec(1) * vec(1)) < Parameters::max_radius2 && fabs(vec(2)) < Parameters::max_z); - } - - template - constexpr float magneticFieldAtPoint(Vec3 const& vec) { - if (isValid(vec)) { - return B0Z(vec) * Kr(vec); - } else { - return 0; - } - } - -} // namespace magneticFieldParabolicPortable - -#endif diff --git a/MagneticField/Portable/BuildFile.xml b/MagneticField/Portable/BuildFile.xml new file mode 100644 index 0000000000000..d186d5ab62fe2 --- /dev/null +++ b/MagneticField/Portable/BuildFile.xml @@ -0,0 +1 @@ + diff --git a/MagneticField/Portable/README.md b/MagneticField/Portable/README.md new file mode 100644 index 0000000000000..6a32bdb37f9cb --- /dev/null +++ b/MagneticField/Portable/README.md @@ -0,0 +1,14 @@ + + # MagneticField/Portable + +This package will host simple and efficient implementations of the magnetic field map for use in well-defined scenarios on GPUs and other heterogeneous hardware. + + ## ParabolicMagneticField + +The parabolic parametrisation of the magnetic field was optimised for extrapolation within R < 115 cm and |Z| < 280 cm, which is (almost) up to the surface of ECAL. +Beyond that it will diverge from the actual magnetic field map, especially at large eta, in ways that are hard to quantify. + +This is based on a simple parabolic formula of Bz(z), neglecting the Br and Bphi components, and is thus vert fast. To account for the missing Br component, which contributes to transverse bending, the parametrization is tuned to reproduce as closely as possible the field integral computed over a straight path originated in (0,0,0) and across the whole tracker, rather than the value of Bz. More details can be found [here](https://twiki.cern.ch/twiki/bin/view/CMSPublic/SWGuideMagneticField#Alternative_parametrizations_of). + + + diff --git a/MagneticField/Portable/interface/ParabolicMagneticField.h b/MagneticField/Portable/interface/ParabolicMagneticField.h new file mode 100644 index 0000000000000..46b88401e58ab --- /dev/null +++ b/MagneticField/Portable/interface/ParabolicMagneticField.h @@ -0,0 +1,44 @@ +#ifndef MagneticField_Portable_interface_ParabolicMagneticField_h +#define MagneticField_Portable_interface_ParabolicMagneticField_h + +#include + +namespace portableParabolicMagneticField { + + struct Parameters { + // These parameters are the best fit of 3.8T to the OAEParametrizedMagneticField parametrization. + + static constexpr float c1 = 3.8114; + static constexpr float b0 = -3.94991e-06; + static constexpr float b1 = 7.53701e-06; + static constexpr float a = 2.43878e-11; + static constexpr float max_radius2 = 13225.f; // tracker radius + static constexpr float max_z = 280.f; // tracker z + }; + + template + constexpr float Kr(Vec3 const& vec) { + return Parameters::a * (vec[0] * vec[0] + vec[1] * vec[1]) + 1.f; + } + + template + constexpr float B0Z(Vec3 const& vec) { + return Parameters::b0 * vec[2] * vec[2] + Parameters::b1 * vec[2] + Parameters::c1; + } + + template + constexpr inline bool isValid(Vec3 const& vec) { + return ((vec[0] * vec[0] + vec[1] * vec[1]) < Parameters::max_radius2 && xtd::abs(vec[2]) < Parameters::max_z); + } + + template + constexpr inline float magneticFieldAtPoint(Vec3 const& vec) { + if (isValid(vec)) { + return B0Z(vec) * Kr(vec); + } else { + return 0; + } + } + +} // namespace portableParabolicMagneticField +#endif diff --git a/MagneticField/ParametrizedEngine/test/BuildFile.xml b/MagneticField/Portable/test/BuildFile.xml similarity index 61% rename from MagneticField/ParametrizedEngine/test/BuildFile.xml rename to MagneticField/Portable/test/BuildFile.xml index a857ea5ca4d74..2c64a855771cb 100644 --- a/MagneticField/ParametrizedEngine/test/BuildFile.xml +++ b/MagneticField/Portable/test/BuildFile.xml @@ -1,9 +1,9 @@ - + - + diff --git a/MagneticField/ParametrizedEngine/test/alpaka/testParabolicParametrizedMagneticField.dev.cc b/MagneticField/Portable/test/alpaka/testParabolicParametrizedTrackerField.dev.cc similarity index 92% rename from MagneticField/ParametrizedEngine/test/alpaka/testParabolicParametrizedMagneticField.dev.cc rename to MagneticField/Portable/test/alpaka/testParabolicParametrizedTrackerField.dev.cc index 8dbb08ac51d9f..0154bfa06b63a 100644 --- a/MagneticField/ParametrizedEngine/test/alpaka/testParabolicParametrizedMagneticField.dev.cc +++ b/MagneticField/Portable/test/alpaka/testParabolicParametrizedTrackerField.dev.cc @@ -11,15 +11,15 @@ #include "HeterogeneousCore/AlpakaInterface/interface/config.h" #include "HeterogeneousCore/AlpakaInterface/interface/memory.h" #include "HeterogeneousCore/AlpakaInterface/interface/workdivision.h" -#include "MagneticField/ParametrizedEngine/interface/ParabolicParametrizedMagneticField.h" +#include "MagneticField/Portable/interface/ParabolicMagneticField.h" using namespace edm; using namespace std; using namespace ALPAKA_ACCELERATOR_NAMESPACE; -using namespace magneticFieldParabolicPortable; +using namespace portableParabolicMagneticField; using Vector3f = Eigen::Matrix; -struct MagneticFieldKernel { +struct TrackerMagneticFieldKernel { template ALPAKA_FN_ACC void operator()(Acc1D const& acc, T const* __restrict__ in, T* __restrict__ out, size_t size) const { for (auto index : cms::alpakatools::uniform_elements(acc, size)) { @@ -77,7 +77,8 @@ int main() { field_host[i] = Vector3f::Zero(); } - float resolution = 0.2; + const float resolution = 0.2; + float maxdelta = 0.; int fail = 0; @@ -94,7 +95,7 @@ int main() { alpaka::memset(queue, field_dev, 0.); auto workDiv = cms::alpakatools::make_workdiv(1, size); - alpaka::exec(queue, workDiv, MagneticFieldKernel{}, points_dev.data(), field_dev.data(), size); + alpaka::exec(queue, workDiv, TrackerMagneticFieldKernel{}, points_dev.data(), field_dev.data(), size); // copy the results from the device to the host alpaka::memcpy(queue, field_host, field_dev); @@ -107,7 +108,7 @@ int main() { const auto& point = points[i]; const auto& referenceB = referenceB_vec[i]; GlobalVector parametricB(field_host[i](0), field_host[i](1), field_host[i](2)); - float delta = (referenceB - parametricB).mag(); + const float delta = (referenceB - parametricB).mag(); if (delta > resolution) { ++fail; if (delta > maxdelta)