Make the magnetic field evaluation reentrant

This removes the shared scratch state from the Chebyshev field evaluation, so that one field can be evaluated from several threads.

- Chebyshev3D keeps the mapped coordinates in a local array instead of a member.
- Chebyshev3DCalc evaluates into scratch on the stack, or into a thread-local buffer for parameterisations larger than 64 rows or columns. The ALICE maps use at most 14.
- MagneticWrapperChebyshev no longer uses static coordinates in the TPC integrals.
- The arithmetic is unchanged, so the field values are bit-identical.
- A unit test evaluates the field from 8 threads and compares with a sequential pass.

This prepares the simulation for multithreaded Geant4, where all workers share one field.

https://its.cern.ch/jira/browse/O2-7305

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

Author: sawenzel

Common/Field/src/MagneticWrapperChebyshev.cxx

@@ -690,7 +690,7 @@ Int_t MagneticWrapperChebyshev::findTPCRatSegment(const Double_t* rpz) const
 
 void MagneticWrapperChebyshev::getTPCIntegral(const Double_t* xyz, Double_t* b) const
 {
-  static Double_t rphiz[3];
+  Double_t rphiz[3];
 
   // TPCInt region
   // convert coordinates to cyl system
@@ -712,7 +712,7 @@ void MagneticWrapperChebyshev::getTPCIntegral(const Double_t* xyz, Double_t* b)
 
 void MagneticWrapperChebyshev::getTPCRatIntegral(const Double_t* xyz, Double_t* b) const
 {
-  static Double_t rphiz[3];
+  Double_t rphiz[3];
 
   // TPCRatIntegral region
   // convert coordinates to cylindrical system

Common/Field/test/testMagneticField.cxx

@@ -21,6 +21,10 @@
 #include <TFile.h>
 #include <TStopwatch.h>
 #include <TRandom.h>
+#include <algorithm>
+#include <array>
+#include <thread>
+#include <vector>
 
 using namespace o2::field;
 
@@ -144,3 +148,40 @@ BOOST_AUTO_TEST_CASE(MagneticField_reinitialization_test)
     }
   }
 }
+
+BOOST_AUTO_TEST_CASE(MagneticField_threads)
+{
+  // evaluating one field from several threads at once gives bit-identical values to a sequential pass
+  MagneticField fld("Maps", "Maps", 1., 1., o2::field::MagFieldParam::k5kG);
+  constexpr int npoints = 20000;
+  constexpr int nthreads = 8;
+  std::vector<std::array<double, 3>> xyz(npoints);
+  for (auto& p : xyz) {
+    p = {gRandom->Uniform(-450., 450.), gRandom->Uniform(-450., 450.), gRandom->Uniform(-1500., 1500.)};
+  }
+  // field, Bz, TPC integral and TPC ratio integral for every point
+  constexpr int nval = 10;
+  auto evaluate = [&fld, &xyz](std::vector<double>& out) {
+    out.resize(npoints * nval);
+    for (int i = 0; i < npoints; ++i) {
+      fld.Field(xyz[i].data(), &out[i * nval]);
+      out[i * nval + 3] = fld.getBz(xyz[i].data());
+      fld.getTPCIntegral(xyz[i].data(), &out[i * nval + 4]);
+      fld.getTPCRatIntegral(xyz[i].data(), &out[i * nval + 7]);
+    }
+  };
+  std::vector<double> reference;
+  evaluate(reference);
+  std::vector<std::vector<double>> results(nthreads);
+  std::vector<std::thread> threads;
+  for (int t = 0; t < nthreads; ++t) {
+    threads.emplace_back(evaluate, std::ref(results[t]));
+  }
+  for (auto& th : threads) {
+    th.join();
+  }
+  for (int t = 0; t < nthreads; ++t) {
+    auto diff = std::mismatch(reference.begin(), reference.end(), results[t].begin());
+    BOOST_CHECK_MESSAGE(diff.first == reference.end(), "thread " << t << " differs at value " << (diff.first - reference.begin()));
+  }
+}

Common/MathUtils/include/MathUtils/Chebyshev3D.h

@@ -292,66 +292,72 @@ inline Bool_t Chebyshev3D::isInside(const Double_t* par) const
 /// Evaluates Chebyshev parameterization for 3d->DimOut function
 inline void Chebyshev3D::Eval(const Float_t* par, Float_t* res)
 {
+  Float_t coef[3]; // mapped arguments, local so that evaluation is reentrant
   for (int i = 3; i--;) {
-    mTemporaryCoefficient[i] = mapToInternal(par[i], i);
+    coef[i] = mapToInternal(par[i], i);
   }
   for (int i = mOutputArrayDimension; i--;) {
-    res[i] = getChebyshevCalc(i)->Eval(mTemporaryCoefficient);
+    res[i] = getChebyshevCalc(i)->Eval(coef);
   }
 }
 
 /// Evaluates Chebyshev parameterization for 3d->DimOut function
 inline void Chebyshev3D::Eval(const Double_t* par, Double_t* res)
 {
+  Float_t coef[3]; // mapped arguments, local so that evaluation is reentrant
   for (int i = 3; i--;) {
-    mTemporaryCoefficient[i] = mapToInternal(par[i], i);
+    coef[i] = mapToInternal(par[i], i);
   }
   for (int i = mOutputArrayDimension; i--;) {
-    res[i] = getChebyshevCalc(i)->Eval(mTemporaryCoefficient);
+    res[i] = getChebyshevCalc(i)->Eval(coef);
   }
 }
 
 /// Evaluates Chebyshev parameterization for idim-th output dimension of 3d->DimOut function
 inline Double_t Chebyshev3D::Eval(const Double_t* par, int idim)
 {
+  Float_t coef[3]; // mapped arguments, local so that evaluation is reentrant
   for (int i = 3; i--;) {
-    mTemporaryCoefficient[i] = mapToInternal(par[i], i);
+    coef[i] = mapToInternal(par[i], i);
   }
-  return getChebyshevCalc(idim)->Eval(mTemporaryCoefficient);
+  return getChebyshevCalc(idim)->Eval(coef);
 }
 
 /// Evaluates Chebyshev parameterization for idim-th output dimension of 3d->DimOut function
 inline Float_t Chebyshev3D::Eval(const Float_t* par, int idim)
 {
+  Float_t coef[3]; // mapped arguments, local so that evaluation is reentrant
   for (int i = 3; i--;) {
-    mTemporaryCoefficient[i] = mapToInternal(par[i], i);
+    coef[i] = mapToInternal(par[i], i);
   }
-  return getChebyshevCalc(idim)->Eval(mTemporaryCoefficient);
+  return getChebyshevCalc(idim)->Eval(coef);
 }
 
 /// Returns the gradient matrix
 inline void Chebyshev3D::evaluateDerivative3D(const Float_t* par, Float_t dbdr[3][3])
 {
+  Float_t coef[3]; // mapped arguments, local so that evaluation is reentrant
   for (int i = 3; i--;) {
-    mTemporaryCoefficient[i] = mapToInternal(par[i], i);
+    coef[i] = mapToInternal(par[i], i);
   }
   for (int ib = 3; ib--;) {
     for (int id = 3; id--;) {
-      dbdr[ib][id] = getChebyshevCalc(ib)->evaluateDerivative(id, mTemporaryCoefficient) * mBoundaryMappingScale[id];
+      dbdr[ib][id] = getChebyshevCalc(ib)->evaluateDerivative(id, coef) * mBoundaryMappingScale[id];
     }
   }
 }
 
 /// Returns the gradient matrix
 inline void Chebyshev3D::evaluateDerivative3D2(const Float_t* par, Float_t dbdrdr[3][3][3])
 {
+  Float_t coef[3]; // mapped arguments, local so that evaluation is reentrant
   for (int i = 3; i--;) {
-    mTemporaryCoefficient[i] = mapToInternal(par[i], i);
+    coef[i] = mapToInternal(par[i], i);
   }
   for (int ib = 3; ib--;) {
     for (int id = 3; id--;) {
       for (int id1 = 3; id1--;) {
-        dbdrdr[ib][id][id1] = getChebyshevCalc(ib)->evaluateDerivative2(id, id1, mTemporaryCoefficient) *
+        dbdrdr[ib][id][id1] = getChebyshevCalc(ib)->evaluateDerivative2(id, id1, coef) *
                               mBoundaryMappingScale[id] * mBoundaryMappingScale[id1];
       }
     }
@@ -361,22 +367,24 @@ inline void Chebyshev3D::evaluateDerivative3D2(const Float_t* par, Float_t dbdrd
 // Evaluates Chebyshev parameterization derivative for 3d->DimOut function
 inline void Chebyshev3D::evaluateDerivative(int dimd, const Float_t* par, Float_t* res)
 {
+  Float_t coef[3]; // mapped arguments, local so that evaluation is reentrant
   for (int i = 3; i--;) {
-    mTemporaryCoefficient[i] = mapToInternal(par[i], i);
+    coef[i] = mapToInternal(par[i], i);
   }
   for (int i = mOutputArrayDimension; i--;) {
-    res[i] = getChebyshevCalc(i)->evaluateDerivative(dimd, mTemporaryCoefficient) * mBoundaryMappingScale[dimd];
+    res[i] = getChebyshevCalc(i)->evaluateDerivative(dimd, coef) * mBoundaryMappingScale[dimd];
   };
 }
 
 // Evaluates Chebyshev parameterization 2nd derivative over dimd1 and dimd2 dimensions for 3d->DimOut function
 inline void Chebyshev3D::evaluateDerivative2(int dimd1, int dimd2, const Float_t* par, Float_t* res)
 {
+  Float_t coef[3]; // mapped arguments, local so that evaluation is reentrant
   for (int i = 3; i--;) {
-    mTemporaryCoefficient[i] = mapToInternal(par[i], i);
+    coef[i] = mapToInternal(par[i], i);
   }
   for (int i = mOutputArrayDimension; i--;) {
-    res[i] = getChebyshevCalc(i)->evaluateDerivative2(dimd1, dimd2, mTemporaryCoefficient) *
+    res[i] = getChebyshevCalc(i)->evaluateDerivative2(dimd1, dimd2, coef) *
              mBoundaryMappingScale[dimd1] * mBoundaryMappingScale[dimd2];
   }
 }
@@ -385,20 +393,22 @@ inline void Chebyshev3D::evaluateDerivative2(int dimd1, int dimd2, const Float_t
 /// function
 inline Float_t Chebyshev3D::evaluateDerivative(int dimd, const Float_t* par, int idim)
 {
+  Float_t coef[3]; // mapped arguments, local so that evaluation is reentrant
   for (int i = 3; i--;) {
-    mTemporaryCoefficient[i] = mapToInternal(par[i], i);
+    coef[i] = mapToInternal(par[i], i);
   }
-  return getChebyshevCalc(idim)->evaluateDerivative(dimd, mTemporaryCoefficient) * mBoundaryMappingScale[dimd];
+  return getChebyshevCalc(idim)->evaluateDerivative(dimd, coef) * mBoundaryMappingScale[dimd];
 }
 
 /// Evaluates Chebyshev parameterization 2ns derivative over dimd1 and dimd2 dimensions for idim-th output dimension of
 /// 3d->DimOut function
 inline Float_t Chebyshev3D::evaluateDerivative2(int dimd1, int dimd2, const Float_t* par, int idim)
 {
+  Float_t coef[3]; // mapped arguments, local so that evaluation is reentrant
   for (int i = 3; i--;) {
-    mTemporaryCoefficient[i] = mapToInternal(par[i], i);
+    coef[i] = mapToInternal(par[i], i);
   }
-  return getChebyshevCalc(idim)->evaluateDerivative2(dimd1, dimd2, mTemporaryCoefficient) *
+  return getChebyshevCalc(idim)->evaluateDerivative2(dimd1, dimd2, coef) *
          mBoundaryMappingScale[dimd1] * mBoundaryMappingScale[dimd2];
 }
 

Common/MathUtils/include/MathUtils/Chebyshev3DCalc.h

@@ -19,6 +19,8 @@
 #include <TNamed.h> // for TNamed
 #include <cstdio>   // for FILE, stdout
 #include <cmath>    // for std::fma
+#include <algorithm>
+#include <vector>
 #include "Rtypes.h" // for Float_t, UShort_t, Int_t, Double_t, etc
 
 class TString;
@@ -171,6 +173,37 @@ class Chebyshev3DCalc : public TNamed
   Double_t Eval(const Double_t* par) const;
 
  private:
+  /// Scratch arrays for one evaluation: on the stack for the usual sizes, thread-local beyond them,
+  /// so that evaluation is reentrant. The ALICE field maps use at most 14 rows and columns.
+  class EvalScratch
+  {
+   public:
+    EvalScratch(int nrows, int ncols)
+    {
+      if (nrows <= kStackSize && ncols <= kStackSize) {
+        mRows = mStackRows;
+        mCols = mStackCols;
+      } else {
+        thread_local std::vector<Float_t> rows, cols;
+        rows.resize(std::max<size_t>(nrows, rows.size()));
+        cols.resize(std::max<size_t>(ncols, cols.size()));
+        mRows = rows.data();
+        mCols = cols.data();
+      }
+    }
+    EvalScratch(const EvalScratch&) = delete;
+    EvalScratch& operator=(const EvalScratch&) = delete;
+    Float_t* rows() { return mRows; }
+    Float_t* cols() { return mCols; }
+
+   private:
+    static constexpr int kStackSize = 64;
+    Float_t mStackRows[kStackSize];
+    Float_t mStackCols[kStackSize];
+    Float_t* mRows;
+    Float_t* mCols;
+  };
+
   Int_t mNumberOfCoefficients;    ///< total number of coeeficients
   Int_t mNumberOfRows;            ///< number of significant rows in the 3D coeffs matrix
   Int_t mNumberOfColumns;         ///< max number of significant cols in the 3D coeffs matrix
@@ -188,8 +221,8 @@ class Chebyshev3DCalc : public TNamed
   // coeffs for col/row
   Float_t* mCoefficients; //[mNumberOfCoefficients] array of Chebyshev coefficients
 
-  Float_t* mTemporaryCoefficients2D; //[mNumberOfColumns] temp. coeffs for 2d summation
-  Float_t* mTemporaryCoefficients1D; //[mNumberOfRows] temp. coeffs for 1d summation
+  Float_t* mTemporaryCoefficients2D; //[mNumberOfColumns] no longer used for evaluation
+  Float_t* mTemporaryCoefficients1D; //[mNumberOfRows] no longer used for evaluation
 
   ClassDefOverride(o2::math_utils::Chebyshev3DCalc,
                    2) // Class for interpolation of 3D->1 function by Chebyshev parametrization
@@ -223,32 +256,38 @@ inline Float_t Chebyshev3DCalc::chebyshevEvaluation1D(Float_t x, const Float_t*
 /// VERY IMPORTANT: par must contain the function arguments ALREADY MAPPED to [-1:1] interval
 inline Float_t Chebyshev3DCalc::Eval(const Float_t* par) const
 {
+  EvalScratch scratch(mNumberOfRows, mNumberOfColumns);
+  Float_t* coefs2D = scratch.cols();
+  Float_t* coefs1D = scratch.rows();
   for (int id0 = mNumberOfRows; id0--;) {
     int nCLoc = mNumberOfColumnsAtRow[id0]; // number of significant coefs on this row
     int col0 = mColumnAtRowBeginning[id0];  // beginning of local column in the 2D boundary matrix
     for (int id1 = nCLoc; id1--;) {
       int id = id1 + col0;
-      mTemporaryCoefficients2D[id1] = chebyshevEvaluation1D(par[2], mCoefficients + mCoefficientBound2D1[id], mCoefficientBound2D0[id]);
+      coefs2D[id1] = chebyshevEvaluation1D(par[2], mCoefficients + mCoefficientBound2D1[id], mCoefficientBound2D0[id]);
     }
-    mTemporaryCoefficients1D[id0] = chebyshevEvaluation1D(par[1], mTemporaryCoefficients2D, nCLoc);
+    coefs1D[id0] = chebyshevEvaluation1D(par[1], coefs2D, nCLoc);
   }
-  return chebyshevEvaluation1D(par[0], mTemporaryCoefficients1D, mNumberOfRows);
+  return chebyshevEvaluation1D(par[0], coefs1D, mNumberOfRows);
 }
 
 /// Evaluates Chebyshev parameterization for 3D function.
 /// VERY IMPORTANT: par must contain the function arguments ALREADY MAPPED to [-1:1] interval
 inline Double_t Chebyshev3DCalc::Eval(const Double_t* par) const
 {
+  EvalScratch scratch(mNumberOfRows, mNumberOfColumns);
+  Float_t* coefs2D = scratch.cols();
+  Float_t* coefs1D = scratch.rows();
   for (int id0 = mNumberOfRows; id0--;) {
     int nCLoc = mNumberOfColumnsAtRow[id0]; // number of significant coefs on this row
     int col0 = mColumnAtRowBeginning[id0];  // beginning of local column in the 2D boundary matrix
     for (int id1 = nCLoc; id1--;) {
       int id = id1 + col0;
-      mTemporaryCoefficients2D[id1] = chebyshevEvaluation1D(par[2], mCoefficients + mCoefficientBound2D1[id], mCoefficientBound2D0[id]);
+      coefs2D[id1] = chebyshevEvaluation1D(par[2], mCoefficients + mCoefficientBound2D1[id], mCoefficientBound2D0[id]);
     }
-    mTemporaryCoefficients1D[id0] = chebyshevEvaluation1D(par[1], mTemporaryCoefficients2D, nCLoc);
+    coefs1D[id0] = chebyshevEvaluation1D(par[1], coefs2D, nCLoc);
   }
-  return chebyshevEvaluation1D(par[0], mTemporaryCoefficients1D, mNumberOfRows);
+  return chebyshevEvaluation1D(par[0], coefs1D, mNumberOfRows);
 }
 } // namespace math_utils
 } // namespace o2