diff --git a/Common/Field/include/Field/MagneticField.h b/Common/Field/include/Field/MagneticField.h index 7c95feacb8787..14a0bff30ae08 100644 --- a/Common/Field/include/Field/MagneticField.h +++ b/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 diff --git a/Common/Field/src/MagneticField.cxx b/Common/Field/src/MagneticField.cxx index 768e42bef6875..603a89d160ba0 100644 --- a/Common/Field/src/MagneticField.cxx +++ b/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"); diff --git a/Common/Field/test/testMagneticField.cxx b/Common/Field/test/testMagneticField.cxx index 9fa8c92260458..df300a4ac7194 100644 --- a/Common/Field/test/testMagneticField.cxx +++ b/Common/Field/test/testMagneticField.cxx @@ -18,6 +18,7 @@ #include "Field/MagFieldFast.h" #include #include // for FairLogger +#include #include #include @@ -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 fld = std::make_unique("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("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); + } + } +}