Alpaka magnetic field polish

MagneticField/Portable/BuildFile.xml

@@ -0,0 +1 @@
+<use name="xtd"/>

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).
+
+
+

MagneticField/Portable/interface/ParabolicMagneticField.h

@@ -0,0 +1,44 @@
+#ifndef MagneticField_Portable_interface_ParabolicMagneticField_h
+#define MagneticField_Portable_interface_ParabolicMagneticField_h
+
+#include <xtd/stdlib/abs.h>
+
+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 <typename Vec3>
+  constexpr float Kr(Vec3 const& vec) {
+    return Parameters::a * (vec[0] * vec[0] + vec[1] * vec[1]) + 1.f;
+  }
+
+  template <typename Vec3>
+  constexpr float B0Z(Vec3 const& vec) {
+    return Parameters::b0 * vec[2] * vec[2] + Parameters::b1 * vec[2] + Parameters::c1;
+  }
+
+  template <typename Vec3>
+  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 <typename Vec3>
+  constexpr inline float parabolicMagneticFieldAtPoint(Vec3 const& vec) {
+    if (isValid(vec)) {
+      return B0Z(vec) * Kr(vec);
+    } else {
+      return 0;
+    }
+  }
+
+}  // namespace portableParabolicMagneticField
+#endif

MagneticField/ParametrizedEngine/test/BuildFile.xml → MagneticField/Portable/test/BuildFile.xml

@@ -1,9 +1,9 @@
-<bin file="alpaka/testParabolicParametrizedMagneticField.dev.cc">
+<bin file="alpaka/testParabolicParametrizedTrackerField.dev.cc">
   <use name="alpaka"/>
   <use name="eigen"/>
   <use name="DataFormats/GeometryVector"/>
   <use name="FWCore/Utilities"/>
   <use name="HeterogeneousCore/AlpakaInterface"/>
-  <use name="MagneticField/ParametrizedEngine" source_only="1"/>
+  <use name="MagneticField/Portable"/>
   <flags ALPAKA_BACKENDS="1"/>
 </bin>

MagneticField/ParametrizedEngine/test/alpaka/testParabolicParametrizedMagneticField.dev.cc → 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<float, 3, 1>;
 
-struct MagneticFieldKernel {
+struct TrackerMagneticFieldKernel {
   template <typename T>
   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<Acc1D>(1, size);
-    alpaka::exec<Acc1D>(queue, workDiv, MagneticFieldKernel{}, points_dev.data(), field_dev.data(), size);
+    alpaka::exec<Acc1D>(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)