From e0c02526f7c6cd62c0339dff8e7c84b65bb26b01 Mon Sep 17 00:00:00 2001 From: alexstrel Date: Tue, 12 Aug 2025 13:51:34 -0500 Subject: [PATCH 1/4] rebased to the latest IB --- FWCore/Utilities/interface/abs.h | 22 +++++++++ .../ParabolicParametrizedMagneticField.h | 47 ------------------- MagneticField/Portable/README.md | 14 ++++++ .../interface/ParabolicMagneticField.h | 44 +++++++++++++++++ .../test/BuildFile.xml | 4 +- ...tParabolicParametrizedTrackerField.dev.cc} | 15 +++--- 6 files changed, 90 insertions(+), 56 deletions(-) create mode 100644 FWCore/Utilities/interface/abs.h delete mode 100644 MagneticField/ParametrizedEngine/interface/ParabolicParametrizedMagneticField.h create mode 100644 MagneticField/Portable/README.md create mode 100644 MagneticField/Portable/interface/ParabolicMagneticField.h rename MagneticField/{ParametrizedEngine => Portable}/test/BuildFile.xml (61%) rename MagneticField/{ParametrizedEngine/test/alpaka/testParabolicParametrizedMagneticField.dev.cc => Portable/test/alpaka/testParabolicParametrizedTrackerField.dev.cc} (91%) diff --git a/FWCore/Utilities/interface/abs.h b/FWCore/Utilities/interface/abs.h new file mode 100644 index 0000000000000..39fe1fb2a965c --- /dev/null +++ b/FWCore/Utilities/interface/abs.h @@ -0,0 +1,22 @@ +#ifndef FWCore_Utilities_Interface_abs_h +#define FWCore_Utilities_Interface_abs_h + +#include + +namespace edm { + + template + concept arithmetic = std::integral || std::floating_point; + + template + constexpr T abs(T x) noexcept { + if constexpr (std::unsigned_integral) { + return x; + } else { + return x < T{0} ? T{-x} : x; + } + } + +} // namespace edm + +#endif //FWCore_Utilities_Interface_abs_h 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/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..9a2806345689b --- /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 && edm::abs(acc, vec[2]) < Parameters::max_z); + } + + template + constexpr inline float parabolicMagneticFieldAtPoint(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..607c98c202a3c 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 91% rename from MagneticField/ParametrizedEngine/test/alpaka/testParabolicParametrizedMagneticField.dev.cc rename to MagneticField/Portable/test/alpaka/testParabolicParametrizedTrackerField.dev.cc index 8dbb08ac51d9f..f630375228be9 100644 --- a/MagneticField/ParametrizedEngine/test/alpaka/testParabolicParametrizedMagneticField.dev.cc +++ b/MagneticField/Portable/test/alpaka/testParabolicParametrizedTrackerField.dev.cc @@ -11,21 +11,21 @@ #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)) { out[index](0) = 0; out[index](1) = 0; - out[index](2) = magneticFieldAtPoint(in[index]); + out[index](2) = parabolicMagneticFieldAtPoint(in[index]); } } }; @@ -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) From 49340564cdc1901b7064559a4cd1428eeb7b5e24 Mon Sep 17 00:00:00 2001 From: Alexei Strelchenko Date: Mon, 10 Nov 2025 09:46:32 -0600 Subject: [PATCH 2/4] introduced (replaced) abs with xtd::abs --- FWCore/Utilities/interface/abs.h | 22 ------------------- MagneticField/Portable/BuildFile.xml | 1 + .../interface/ParabolicMagneticField.h | 4 ++-- MagneticField/Portable/test/BuildFile.xml | 1 - 4 files changed, 3 insertions(+), 25 deletions(-) delete mode 100644 FWCore/Utilities/interface/abs.h create mode 100644 MagneticField/Portable/BuildFile.xml diff --git a/FWCore/Utilities/interface/abs.h b/FWCore/Utilities/interface/abs.h deleted file mode 100644 index 39fe1fb2a965c..0000000000000 --- a/FWCore/Utilities/interface/abs.h +++ /dev/null @@ -1,22 +0,0 @@ -#ifndef FWCore_Utilities_Interface_abs_h -#define FWCore_Utilities_Interface_abs_h - -#include - -namespace edm { - - template - concept arithmetic = std::integral || std::floating_point; - - template - constexpr T abs(T x) noexcept { - if constexpr (std::unsigned_integral) { - return x; - } else { - return x < T{0} ? T{-x} : x; - } - } - -} // namespace edm - -#endif //FWCore_Utilities_Interface_abs_h 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/interface/ParabolicMagneticField.h b/MagneticField/Portable/interface/ParabolicMagneticField.h index 9a2806345689b..0ce3f624c71ed 100644 --- a/MagneticField/Portable/interface/ParabolicMagneticField.h +++ b/MagneticField/Portable/interface/ParabolicMagneticField.h @@ -1,7 +1,7 @@ #ifndef MagneticField_Portable_interface_ParabolicMagneticField_h #define MagneticField_Portable_interface_ParabolicMagneticField_h -#include +#include namespace portableParabolicMagneticField { @@ -28,7 +28,7 @@ namespace portableParabolicMagneticField { template constexpr inline bool isValid(Vec3 const& vec) { - return ((vec[0] * vec[0] + vec[1] * vec[1]) < Parameters::max_radius2 && edm::abs(acc, vec[2]) < Parameters::max_z); + return ((vec[0] * vec[0] + vec[1] * vec[1]) < Parameters::max_radius2 && xtd::abs(vec[2]) < Parameters::max_z); } template diff --git a/MagneticField/Portable/test/BuildFile.xml b/MagneticField/Portable/test/BuildFile.xml index 607c98c202a3c..55437b02470f2 100644 --- a/MagneticField/Portable/test/BuildFile.xml +++ b/MagneticField/Portable/test/BuildFile.xml @@ -4,6 +4,5 @@ - From 96c0e68d14df449b84df054b2b1eaab028717008 Mon Sep 17 00:00:00 2001 From: Alexei Strelchenko Date: Mon, 10 Nov 2025 11:00:04 -0600 Subject: [PATCH 3/4] corrected BuildFile.xml --- MagneticField/Portable/test/BuildFile.xml | 1 + 1 file changed, 1 insertion(+) diff --git a/MagneticField/Portable/test/BuildFile.xml b/MagneticField/Portable/test/BuildFile.xml index 55437b02470f2..2c64a855771cb 100644 --- a/MagneticField/Portable/test/BuildFile.xml +++ b/MagneticField/Portable/test/BuildFile.xml @@ -4,5 +4,6 @@ + From 739777161eab9576b958c9916e3553e8bbde0575 Mon Sep 17 00:00:00 2001 From: Alexei Strelchenko Date: Thu, 13 Nov 2025 09:20:34 -0600 Subject: [PATCH 4/4] restored the original name of a function --- MagneticField/Portable/interface/ParabolicMagneticField.h | 2 +- .../test/alpaka/testParabolicParametrizedTrackerField.dev.cc | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/MagneticField/Portable/interface/ParabolicMagneticField.h b/MagneticField/Portable/interface/ParabolicMagneticField.h index 0ce3f624c71ed..46b88401e58ab 100644 --- a/MagneticField/Portable/interface/ParabolicMagneticField.h +++ b/MagneticField/Portable/interface/ParabolicMagneticField.h @@ -32,7 +32,7 @@ namespace portableParabolicMagneticField { } template - constexpr inline float parabolicMagneticFieldAtPoint(Vec3 const& vec) { + constexpr inline float magneticFieldAtPoint(Vec3 const& vec) { if (isValid(vec)) { return B0Z(vec) * Kr(vec); } else { diff --git a/MagneticField/Portable/test/alpaka/testParabolicParametrizedTrackerField.dev.cc b/MagneticField/Portable/test/alpaka/testParabolicParametrizedTrackerField.dev.cc index f630375228be9..0154bfa06b63a 100644 --- a/MagneticField/Portable/test/alpaka/testParabolicParametrizedTrackerField.dev.cc +++ b/MagneticField/Portable/test/alpaka/testParabolicParametrizedTrackerField.dev.cc @@ -25,7 +25,7 @@ struct TrackerMagneticFieldKernel { for (auto index : cms::alpakatools::uniform_elements(acc, size)) { out[index](0) = 0; out[index](1) = 0; - out[index](2) = parabolicMagneticFieldAtPoint(in[index]); + out[index](2) = magneticFieldAtPoint(in[index]); } } };