Field: keep polarity on MagneticField re-initialisation

Common/Field/include/Field/MagneticField.h

@@ -262,8 +262,8 @@ class MagneticField : public FairField
 
   Int_t mDefaultIntegration;             ///< Default integration method as indicated in Geant
   Int_t mPrecisionInteg;                 ///< Alternative integration method, e.g. for higher precision
-  Double_t mMultipicativeFactorSolenoid; ///< Multiplicative factor for solenoid
-  Double_t mMultipicativeFactorDipole;   ///< Multiplicative factor for dipole
+  Double_t mMultipicativeFactorSolenoid; ///< Multiplicative factor for solenoid, polarity convention applied
+  Double_t mMultipicativeFactorDipole;   ///< Multiplicative factor for dipole, polarity convention applied
   Double_t mMaxField;                    ///< Max Field as indicated in Geant
   Bool_t mDipoleOnOffFlag;               ///< Dipole ON/OFF flag
 

Common/Field/src/MagneticField.cxx

@@ -108,8 +108,8 @@ MagneticField::MagneticField(const char* name, const char* title, Double_t facto
     mBeamEnergy(be),
     mDefaultIntegration(integ),
     mPrecisionInteg(1),
-    mMultipicativeFactorSolenoid(factorSol),
-    mMultipicativeFactorDipole(factorDip),
+    mMultipicativeFactorSolenoid(1.),
+    mMultipicativeFactorDipole(1.),
     mMaxField(fmax),
     mDipoleOnOffFlag(factorDip == 0.),
     mQuadrupoleGradient(0),
@@ -122,6 +122,8 @@ MagneticField::MagneticField(const char* name, const char* title, Double_t facto
   /*
    * Constructor for human readable params
    */
+  setFactorSolenoid(factorSol);
+  setFactorDipole(factorDip);
   setDataFileName(path.c_str());
   if (!gOriginBias) {
     checkOriginBias();
@@ -139,8 +141,8 @@ MagneticField::MagneticField(const MagFieldParam& param)
     mBeamEnergy(param.GetBeamEnergy()),
     mDefaultIntegration(param.GetDefInt()),
     mPrecisionInteg(1),
-    mMultipicativeFactorSolenoid(param.GetFactorSol()), // temporary
-    mMultipicativeFactorDipole(param.GetFactorDip()),   // temporary
+    mMultipicativeFactorSolenoid(1.),
+    mMultipicativeFactorDipole(1.),
     mMaxField(param.GetMaxField()),
     mDipoleOnOffFlag(param.GetFactorDip() == 0.),
     mQuadrupoleGradient(0),
@@ -153,6 +155,8 @@ MagneticField::MagneticField(const MagFieldParam& param)
   /*
    * Constructor for FairParam derived params
    */
+  setFactorSolenoid(param.GetFactorSol());
+  setFactorDipole(param.GetFactorDip());
   setDataFileName(param.GetMapPath());
   if (!gOriginBias) {
     checkOriginBias();
@@ -235,8 +239,11 @@ void MagneticField::CreateField()
 
   loadParameterization();
   initializeMachineField(mBeamType, mBeamEnergy);
-  setFactorSolenoid(mMultipicativeFactorSolenoid);
-  setFactorDipole(mMultipicativeFactorDipole);
+  // The scaling factors are left alone: they already carry the polarity convention, and
+  // re-applying it would invert a field re-initialized after being read back from a file.
+  if (mFastField) {
+    mFastField->setFactorSol(getFactorSolenoid());
+  }
   double xyz[3] = {0., 0., 0.};
   mSolenoid = getBz(xyz);
   Print("a");

Common/Field/test/testMagneticField.cxx

@@ -18,6 +18,7 @@
 #include "Field/MagFieldFast.h"
 #include <memory>
 #include <fairlogger/Logger.h> // for FairLogger
+#include <TFile.h>
 #include <TStopwatch.h>
 #include <TRandom.h>
 
@@ -98,3 +99,48 @@ BOOST_AUTO_TEST_CASE(MagneticField_test)
     BOOST_CHECK(TMath::Abs(rms[i] / nomBz) < 1.e-3);
   }
 }
+
+BOOST_AUTO_TEST_CASE(MagneticField_reinitialization_test)
+{
+  // The measured map is transient, so a MagneticField read back from a file has to be
+  // re-created before it can be used. That must reproduce the field vectors, not merely
+  // their magnitude: a sign flip leaves |B| untouched.
+  const double points[][3] = {
+    {0., 0., 0.},      // solenoid, on axis
+    {100., 50., 100.}, // solenoid, off axis
+    {10., 10., -900.}, // muon dipole
+    {0., 0., 1000.},   // compensator 1A, side A
+    {0., 0., -2049.},  // compensator 2C, side C
+    {0., 0., 2049.}    // compensator 2A, side A
+  };
+  const int npoints = sizeof(points) / sizeof(points[0]);
+  const double tolerance = 1.e-9; // kGauss
+
+  std::unique_ptr<MagneticField> fld = std::make_unique<MagneticField>("Maps", "Maps", 1., 1., MagFieldParam::k5kG);
+  const double facSol = fld->getFactorSolenoid(), facDip = fld->getFactorDipole();
+  double bref[npoints][3] = {};
+  for (int ip = 0; ip < npoints; ip++) {
+    fld->Field(points[ip], bref[ip]);
+    // a point where the field vanishes would make the comparisons below vacuous
+    BOOST_CHECK(TMath::Abs(bref[ip][0]) + TMath::Abs(bref[ip][1]) + TMath::Abs(bref[ip][2]) > tolerance);
+  }
+
+  const char* fname = "testMagneticFieldReinitialization.root";
+  {
+    TFile fout(fname, "recreate");
+    fout.WriteObject(fld.get(), "field");
+  }
+  TFile fin(fname);
+  auto* fldRead = fin.Get<MagneticField>("field");
+  BOOST_REQUIRE(fldRead != nullptr);
+  fldRead->CreateField();
+  BOOST_CHECK_EQUAL(fldRead->getFactorSolenoid(), facSol);
+  BOOST_CHECK_EQUAL(fldRead->getFactorDipole(), facDip);
+  for (int ip = 0; ip < npoints; ip++) {
+    double b[3] = {};
+    fldRead->Field(points[ip], b);
+    for (int i = 0; i < 3; i++) {
+      BOOST_CHECK_SMALL(b[i] - bref[ip][i], tolerance);
+    }
+  }
+}