diff --git a/Analysis/HitAnalyzer/plugins/ResidualGlobalCorrectionMakerBase.cc b/Analysis/HitAnalyzer/plugins/ResidualGlobalCorrectionMakerBase.cc
index e7c902c8850e3..9e8be2d677f87 100644
--- a/Analysis/HitAnalyzer/plugins/ResidualGlobalCorrectionMakerBase.cc
+++ b/Analysis/HitAnalyzer/plugins/ResidualGlobalCorrectionMakerBase.cc
@@ -76,6 +76,7 @@
 #include "CommonTools/UtilAlgos/interface/TFileService.h"
 
 #include "Alignment/CommonAlignment/interface/Utilities.h"
+#include "TRandom.h"
 
 // #include "DataFormats/Math/interface/Point3D.h"
 
@@ -136,6 +137,8 @@ ResidualGlobalCorrectionMakerBase::ResidualGlobalCorrectionMakerBase(const edm::
   doRes_ = iConfig.getParameter<bool>("doRes");
   useIdealGeometry_ = iConfig.getParameter<bool>("useIdealGeometry");
   corFiles_ = iConfig.getParameter<std::vector<std::string>>("corFiles");
+  fieldlabel_ = iConfig.getParameter<std::string>("MagneticFieldLabel");
+
 
   inputBs_ = consumes<reco::BeamSpot>(edm::InputTag("offlineBeamSpot"));
 
@@ -347,8 +350,63 @@ ResidualGlobalCorrectionMakerBase::beginRun(edm::Run const& run, edm::EventSetup
 //   const MagneticField* field = thePropagator->magneticField();
 
   edm::ESHandle<MagneticField> magfield;
-  es.get<IdealMagneticFieldRecord>().get(magfield);
+  es.get<IdealMagneticFieldRecord>().get(fieldlabel_, magfield);
   auto field = magfield.product();
+
+  constexpr bool dofieldtest = false;
+
+  if (dofieldtest) {
+    constexpr int npoints = 1000;
+
+    std::vector<GlobalPoint> points;
+    points.reserve(npoints);
+
+    std::vector<GlobalVector> fieldfwd(npoints);
+    std::vector<GlobalVector> fieldrev(npoints);
+
+    for (int ipoint = 0; ipoint < npoints; ++ipoint) {
+      // const double x = gRandom->Uniform(-100., 100.);
+      // const double y = gRandom->Uniform(-10., 10.);
+      // const double z = gRandom->Uniform(-10., 10.);
+
+      const double rho = gRandom->Uniform(0., 200.);
+      const double phi = gRandom->Uniform(-M_PI, M_PI);
+      const double z = gRandom->Uniform(-500., 500.);
+
+      const double x = rho*std::cos(phi);
+      const double y = rho*std::sin(phi);
+
+      points.emplace_back(x, y, z);
+    }
+
+    for (int ipoint = 0; ipoint < npoints; ++ipoint) {
+      fieldfwd[ipoint] = field->inTeslaUnchecked(points[ipoint]);
+    }
+
+    for (int ipoint = npoints - 1; ipoint >= 0; --ipoint) {
+      fieldrev[ipoint] = field->inTeslaUnchecked(points[ipoint]);
+    }
+
+    bool fieldmatch = true;
+    for (int ipoint = 0; ipoint < npoints; ++ipoint) {
+      auto const &ifieldfwd = fieldfwd[ipoint];
+      auto const &ifieldrev = fieldrev[ipoint];
+      bool imatch = ifieldfwd.x() == ifieldrev.x() && ifieldfwd.y() == ifieldrev.y() && ifieldfwd.z() == ifieldrev.z();
+      if (!imatch) {
+        std::cout << "field mismatch! ipoint = " << ipoint << " fwd: " << ifieldfwd << " rev: " << ifieldrev << std::endl;
+        fieldmatch = false;
+      }
+    }
+
+    std::cout << "fieldmatch = " << fieldmatch << std::endl;
+
+    if (!fieldmatch) {
+      throw std::logic_error("field values depend on evaluation order");
+    }
+
+
+
+  }
   
   std::set<std::pair<int, DetId> > parmset;
   
diff --git a/Analysis/HitAnalyzer/plugins/ResidualGlobalCorrectionMakerBase.h b/Analysis/HitAnalyzer/plugins/ResidualGlobalCorrectionMakerBase.h
index 922a5c444c897..32b784ccfcffd 100644
--- a/Analysis/HitAnalyzer/plugins/ResidualGlobalCorrectionMakerBase.h
+++ b/Analysis/HitAnalyzer/plugins/ResidualGlobalCorrectionMakerBase.h
@@ -262,6 +262,9 @@ class ResidualGlobalCorrectionMakerBase : public edm::stream::EDProducer<>
   
   
   std::vector<std::string> corFiles_;
+
+  std::string fieldlabel_;
+
   
 //   SiStripClusterInfo siStripClusterInfo_;
 
diff --git a/Analysis/HitAnalyzer/plugins/ResidualGlobalCorrectionMakerTwoTrackG4e.cc b/Analysis/HitAnalyzer/plugins/ResidualGlobalCorrectionMakerTwoTrackG4e.cc
index f4838cfcd664d..6b02b78ab6d9b 100644
--- a/Analysis/HitAnalyzer/plugins/ResidualGlobalCorrectionMakerTwoTrackG4e.cc
+++ b/Analysis/HitAnalyzer/plugins/ResidualGlobalCorrectionMakerTwoTrackG4e.cc
@@ -141,12 +141,16 @@ class ResidualGlobalCorrectionMakerTwoTrackG4e : public ResidualGlobalCorrection
   float Muminusgen_pt;
   float Muminusgen_eta;
   float Muminusgen_phi;
+
+  float Muplusgen_dr;
+  float Muminusgen_dr;
   
   std::array<float, 3> Muplus_refParms;
   std::array<float, 3> Muminus_refParms;
   
   std::vector<float> Muplus_jacRef;
   std::vector<float> Muminus_jacRef;
+  std::vector<float> Jpsi_jacMass;
   
   unsigned int Muplus_nhits;
   unsigned int Muplus_nvalid;
@@ -308,6 +312,9 @@ void ResidualGlobalCorrectionMakerTwoTrackG4e::beginStream(edm::StreamID streami
     tree->Branch("Muminusgen_pt", &Muminusgen_pt);
     tree->Branch("Muminusgen_eta", &Muminusgen_eta);
     tree->Branch("Muminusgen_phi", &Muminusgen_phi);
+
+    tree->Branch("Muplusgen_dr", &Muplusgen_dr);
+    tree->Branch("Muminusgen_dr", &Muminusgen_dr);
     
     tree->Branch("Muplus_refParms", Muplus_refParms.data(), "Muplus_refParms[3]/F");
     tree->Branch("Muminus_refParms", Muminus_refParms.data(), "Muminus_refParms[3]/F");
@@ -315,6 +322,7 @@ void ResidualGlobalCorrectionMakerTwoTrackG4e::beginStream(edm::StreamID streami
     if (fillJac_) {
       tree->Branch("Muplus_jacRef", &Muplus_jacRef);
       tree->Branch("Muminus_jacRef", &Muminus_jacRef);
+      tree->Branch("Jpsi_jacMass", &Jpsi_jacMass);
     }
     
     tree->Branch("Muplus_nhits", &Muplus_nhits);
@@ -535,6 +543,7 @@ void ResidualGlobalCorrectionMakerTwoTrackG4e::produce(edm::Event &iEvent, const
     }
     
     const reco::Candidate *mu0gen = nullptr;
+    double drmin0 = 0.1;
     if (doGen_ && !doSim_) {
       for (auto const &genpart : *genPartCollection) {
         if (genpart.status() != 1) {
@@ -544,9 +553,10 @@ void ResidualGlobalCorrectionMakerTwoTrackG4e::produce(edm::Event &iEvent, const
           continue;
         }
         
-        float dR0 = deltaR(genpart, *itrack);
-        if (dR0 < 0.1 && genpart.charge() == itrack->charge()) {
+        const double dR0 = deltaR(genpart, *itrack);
+        if (dR0 < drmin0 && genpart.charge() == itrack->charge()) {
           mu0gen = &genpart;
+          drmin0 = dR0;
         }
       }
     }
@@ -585,6 +595,7 @@ void ResidualGlobalCorrectionMakerTwoTrackG4e::produce(edm::Event &iEvent, const
       
 
       const reco::Candidate *mu1gen = nullptr;
+      double drmin1 = 0.1;
       
       double massconstraintval = massConstraint_;
       if (doGen_ && !doSim_) {
@@ -596,9 +607,10 @@ void ResidualGlobalCorrectionMakerTwoTrackG4e::produce(edm::Event &iEvent, const
             continue;
           }
           
-          float dR1 = deltaR(genpart, *jtrack);
-          if (dR1 < 0.1 && genpart.charge() == jtrack->charge()) {
+          const double dR1 = deltaR(genpart, *jtrack);
+          if (dR1 < drmin1 && genpart.charge() == jtrack->charge()) {
             mu1gen = &genpart;
+            drmin1 = dR1;
           }
         }
         
@@ -870,6 +882,12 @@ void ResidualGlobalCorrectionMakerTwoTrackG4e::produce(edm::Event &iEvent, const
       bool valid = true;
       
       
+      if (false) {
+        const GlobalPoint fieldrefpoint(itrack->vertex().x(), itrack->vertex().y(), itrack->vertex().z());
+        auto const fieldvalref = field->inTesla(fieldrefpoint);
+        std::cout << "refpos: " << fieldrefpoint << " bfield = " << fieldvalref << std::endl;
+      }
+
       const unsigned int nicons = doMassConstraint_ ? 2 : 1;
 //       const unsigned int nicons = doMassConstraint_ ? 3 : 1;
       
@@ -900,6 +918,7 @@ void ResidualGlobalCorrectionMakerTwoTrackG4e::produce(edm::Event &iEvent, const
         
         if (kinTree->isEmpty() || !kinTree->isConsistent()) {
 //           continue;
+          std::cout << "Abort: invalid kinematic fit!\n";
           valid = false;
           break;
         }
@@ -2202,7 +2221,13 @@ void ResidualGlobalCorrectionMakerTwoTrackG4e::produce(edm::Event &iEvent, const
           const reco::Candidate *muplusgen = nullptr;
           const reco::Candidate *muminusgen = nullptr;
           
+          Muplusgen_dr = -99.;
+          Muminusgen_dr = -99.;
+
           if (doGen_) {
+            double drminplus = 0.1;
+            double drminminus = 0.1;
+
             for (auto const &genpart : *genPartCollection) {
               if (genpart.status() != 1) {
                 continue;
@@ -2210,19 +2235,30 @@ void ResidualGlobalCorrectionMakerTwoTrackG4e::produce(edm::Event &iEvent, const
               if (std::abs(genpart.pdgId()) != 13) {
                 continue;
               }
-              
+
 //               float dRplus = deltaR(genpart.phi(), muarr[idxplus].phi(), genpart.eta(), muarr[idxplus].eta());
-              float dRplus = deltaR(genpart, muarr[idxplus]);
-              if (dRplus < 0.1 && genpart.charge() > 0) {
+              const double dRplus = deltaR(genpart, muarr[idxplus]);
+              if (dRplus < drminplus && genpart.charge() > 0) {
                 muplusgen = &genpart;
+                drminplus = dRplus;
               }
-              
+
 //               float dRminus = deltaR(genpart.phi(), muarr[idxminus].phi(), genpart.eta(), muarr[idxminus].eta());
-              float dRminus = deltaR(genpart, muarr[idxminus]);
-              if (dRminus < 0.1 && genpart.charge() < 0) {
+              const double dRminus = deltaR(genpart, muarr[idxminus]);
+              if (dRminus < drminminus && genpart.charge() < 0) {
                 muminusgen = &genpart;
+                drminminus = dRminus;
               }
             }
+
+            if (muplusgen != nullptr) {
+              Muplusgen_dr = drminplus;
+            }
+
+            if (muminusgen != nullptr) {
+              Muminusgen_dr = drminminus;
+            }
+
           }
           
           if (muplusgen != nullptr) {
@@ -2457,6 +2493,12 @@ void ResidualGlobalCorrectionMakerTwoTrackG4e::produce(edm::Event &iEvent, const
           Muminus_jacRef.resize(3*nparsfinal);
           Map<Matrix<float, 3, Dynamic, RowMajor>>(Muminus_jacRef.data(), 3, nparsfinal) = dxdparms.block(0, trackstateidxminus, nparsfinal, 3).transpose().cast<float>();
 
+          const Matrix<double, 1, 6> mjacalt = massJacobianAltD(refftsarr[0], refftsarr[1], mmu);
+
+          Jpsi_jacMass.resize(nparsfinal);
+          Map<Matrix<float, 1, Dynamic, RowMajor>>(Jpsi_jacMass.data(), 1, nparsfinal) = (mjacalt*dxdparms.leftCols<6>().transpose()).cast<float>();
+
+
           if (false) {
             std::cout << "Muplus pt eta phi = " << Muplus_pt << "  " << Muplus_eta << " " << Muplus_phi << std::endl;
             std::cout << "Muminus pt eta phi = " << Muminus_pt << "  " << Muminus_eta << " " << Muminus_phi << std::endl;
diff --git a/CalibTracker/SiPixelESProducers/interface/SiPixelTemplateDBObjectESProducer.h b/CalibTracker/SiPixelESProducers/interface/SiPixelTemplateDBObjectESProducer.h
index f8699bde4ec94..9c934cf206dd5 100644
--- a/CalibTracker/SiPixelESProducers/interface/SiPixelTemplateDBObjectESProducer.h
+++ b/CalibTracker/SiPixelESProducers/interface/SiPixelTemplateDBObjectESProducer.h
@@ -29,5 +29,8 @@ class SiPixelTemplateDBObjectESProducer : public edm::ESProducer  {
   SiPixelTemplateDBObjectESProducer(const edm::ParameterSet& iConfig);
   ~SiPixelTemplateDBObjectESProducer() override;
   std::shared_ptr<const SiPixelTemplateDBObject> produce(const SiPixelTemplateDBObjectESProducerRcd &);
+
+private:
+  std::string fieldlabel_;
 };
 #endif
diff --git a/CalibTracker/SiPixelESProducers/plugins/SiPixelTemplateDBObjectESProducer.cc b/CalibTracker/SiPixelESProducers/plugins/SiPixelTemplateDBObjectESProducer.cc
index 6f1683d8b5c7c..dfee6796d453c 100644
--- a/CalibTracker/SiPixelESProducers/plugins/SiPixelTemplateDBObjectESProducer.cc
+++ b/CalibTracker/SiPixelESProducers/plugins/SiPixelTemplateDBObjectESProducer.cc
@@ -21,6 +21,7 @@
 using namespace edm;
 
 SiPixelTemplateDBObjectESProducer::SiPixelTemplateDBObjectESProducer(const edm::ParameterSet& iConfig) {
+	fieldlabel_ = iConfig.getParameter<std::string>("MagneticFieldLabel");
 	setWhatProduced(this);
 }
 
@@ -34,7 +35,7 @@ SiPixelTemplateDBObjectESProducer::~SiPixelTemplateDBObjectESProducer(){
 std::shared_ptr<const SiPixelTemplateDBObject> SiPixelTemplateDBObjectESProducer::produce(const SiPixelTemplateDBObjectESProducerRcd & iRecord) {
 	
 	ESHandle<MagneticField> magfield;
-	iRecord.getRecord<IdealMagneticFieldRecord>().get(magfield);
+	iRecord.getRecord<IdealMagneticFieldRecord>().get(fieldlabel_, magfield);
 
 	GlobalPoint center(0.0, 0.0, 0.0);
 	float theMagField = magfield.product()->inTesla(center).mag();
diff --git a/CalibTracker/SiPixelESProducers/python/SiPixelTemplateDBObjectESProducer_cfi.py b/CalibTracker/SiPixelESProducers/python/SiPixelTemplateDBObjectESProducer_cfi.py
index 7d59a51a55301..cdb835f740ff8 100644
--- a/CalibTracker/SiPixelESProducers/python/SiPixelTemplateDBObjectESProducer_cfi.py
+++ b/CalibTracker/SiPixelESProducers/python/SiPixelTemplateDBObjectESProducer_cfi.py
@@ -1,4 +1,6 @@
 import FWCore.ParameterSet.Config as cms
 
-siPixelTemplateDBObjectESProducer = cms.ESProducer("SiPixelTemplateDBObjectESProducer")
+siPixelTemplateDBObjectESProducer = cms.ESProducer("SiPixelTemplateDBObjectESProducer",
+                                                      MagneticFieldLabel = cms.string(""),
+                                                   )
  
diff --git a/MagneticField/ParametrizedEngine/BuildFile.xml b/MagneticField/ParametrizedEngine/BuildFile.xml
index 68e00d770edc5..87ad713c50b2f 100644
--- a/MagneticField/ParametrizedEngine/BuildFile.xml
+++ b/MagneticField/ParametrizedEngine/BuildFile.xml
@@ -1,3 +1,4 @@
+<flags CXXFLAGS="-fno-tree-vectorize"/>
 <use   name="DataFormats/GeometryVector"/>
 <!--use   name="FWCore/Framework"/-->
 <use   name="FWCore/ParameterSet"/>
diff --git a/MagneticField/ParametrizedEngine/src/PolyFit3DParametrizedMagneticField.cc b/MagneticField/ParametrizedEngine/src/PolyFit3DParametrizedMagneticField.cc
index 15ba721ff075f..0f8507a2b8d0b 100644
--- a/MagneticField/ParametrizedEngine/src/PolyFit3DParametrizedMagneticField.cc
+++ b/MagneticField/ParametrizedEngine/src/PolyFit3DParametrizedMagneticField.cc
@@ -63,7 +63,7 @@ PolyFit3DParametrizedMagneticField::inTesla(const GlobalPoint& gp) const {
   if (isDefined(gp)) {
     return inTeslaUnchecked(gp);
   } else {
-    edm::LogWarning("MagneticField|FieldOutsideValidity") << " Point " << gp << " is outside the validity region of PolyFit3DParametrizedMagneticField";
+    // edm::LogWarning("MagneticField|FieldOutsideValidity") << " Point " << gp << " is outside the validity region of PolyFit3DParametrizedMagneticField";
     return GlobalVector();
   }
 }
diff --git a/MagneticField/ParametrizedEngine/src/poly2d_base.cc b/MagneticField/ParametrizedEngine/src/poly2d_base.cc
index f124aa64c4f82..3f876d4f6937a 100644
--- a/MagneticField/ParametrizedEngine/src/poly2d_base.cc
+++ b/MagneticField/ParametrizedEngine/src/poly2d_base.cc
@@ -33,8 +33,8 @@ void poly2d_term::Print(std::ostream &out, bool first_term)
 /////////////////////////////////////////////////////////////////////////////////
 
 //_______________________________________________________________________________
-double     poly2d_base::rval   = 0.; //Last values of r and z used
-double     poly2d_base::zval   = 0.;
+double     poly2d_base::rval   = std::numeric_limits<double>::infinity(); //Last values of r and z used
+double     poly2d_base::zval   = std::numeric_limits<double>::infinity();
 
 double   **poly2d_base::rz_pow = nullptr;  //table with calculated r^n*z^m values
 unsigned   poly2d_base::NTab   = 0;  //rz_pow table size
@@ -57,7 +57,7 @@ poly2d_base::~poly2d_base()
          delete [] rz_pow;
       }
       rz_pow = nullptr;
-      rval = zval = 0.;
+      rval = zval = std::numeric_limits<double>::infinity();
       NPwr = 0;
       NTab = 0;
       rz_set = false;
@@ -81,7 +81,7 @@ void poly2d_base::SetTabSize(const unsigned N)
    for (jr = 1, dN = N; jr < N; ++jr, --dN) {
       rz_pow[jr] = rz_pow[jr-1] + dN;
    }
-   rval = zval = 0.;
+   rval = zval = std::numeric_limits<double>::infinity();
    NTab = N;
 }
 
diff --git a/MagneticField/ParametrizedEngine/src/rz_harm_poly.cc b/MagneticField/ParametrizedEngine/src/rz_harm_poly.cc
index 99ba61a65dd0a..7d1751b7f336c 100644
--- a/MagneticField/ParametrizedEngine/src/rz_harm_poly.cc
+++ b/MagneticField/ParametrizedEngine/src/rz_harm_poly.cc
@@ -12,7 +12,7 @@ using namespace magfieldparam;
 
 //_______________________________________________________________________________
 unsigned rz_harm_poly::Cnt    = 0;      //Number of the "rz_harm_poly" objects
-double   rz_harm_poly::phival = -11111.;//Last phi value used
+double   rz_harm_poly::phival = std::numeric_limits<double>::infinity();//Last phi value used
 bool     rz_harm_poly::phi_set = false; //TRUE if phi value is set
 unsigned rz_harm_poly::MaxM   = 0;      //Max. M among "rz_harm_poly" objects
 
@@ -63,7 +63,7 @@ rz_harm_poly::~rz_harm_poly()
       TrigArr = nullptr;
       TASize = 0;
       MaxM = 0;
-      phival = -11111.;
+      phival = std::numeric_limits<double>::infinity();
       phi_set = false;
    }
 }
@@ -115,6 +115,8 @@ void rz_harm_poly::FillTrigArr(const double phi)
    if (!TrigArr) return;
    trig_pair tp(phi);
    TrigArr[1] = tp;
+   // FIXME this loop will be auto-vectorized in slc7_amd64_gcc700 with default
+   // compiler flags and gives incorrect results
    for (unsigned jp = 2; jp <= MaxM; ++jp) TrigArr[jp] = TrigArr[jp-1].Add(tp);
 }
 
diff --git a/PhysicsTools/NanoAOD/python/muons_cff.py b/PhysicsTools/NanoAOD/python/muons_cff.py
index ad39ee53fe58e..d80058d53126e 100644
--- a/PhysicsTools/NanoAOD/python/muons_cff.py
+++ b/PhysicsTools/NanoAOD/python/muons_cff.py
@@ -144,6 +144,7 @@
                                    applyHitQuality = cms.bool(True),
                                    useIdealGeometry = cms.bool(False),
                                    corFiles = cms.vstring(),
+                                   MagneticFieldLabel = cms.string(""),
 )
 
 trackrefitideal = cms.EDProducer('ResidualGlobalCorrectionMakerG4e',
@@ -165,6 +166,7 @@
                                    applyHitQuality = cms.bool(True),
                                    useIdealGeometry = cms.bool(True),
                                    corFiles = cms.vstring(),
+                                   MagneticFieldLabel = cms.string(""),
 )
 
 mergedGlobalIdxs = cms.EDProducer("GlobalIdxProducer",
diff --git a/PhysicsTools/NanoAOD/python/nano_cff.py b/PhysicsTools/NanoAOD/python/nano_cff.py
index 91e03b1cf1f4d..9939cef9af7bb 100644
--- a/PhysicsTools/NanoAOD/python/nano_cff.py
+++ b/PhysicsTools/NanoAOD/python/nano_cff.py
@@ -451,12 +451,18 @@ def nanoAOD_customizeData(process):
                                     ),
                                 )
 
-    # load 3d field map and use it for g4e propagator
+    # load 3d field map and use it for g4e propagator, geant4 internals via geometry producer and a few other places related to the track refit
     from MagneticField.ParametrizedEngine.parametrizedMagneticField_PolyFit3D_cfi import ParametrizedMagneticFieldProducer as PolyFit3DMagneticFieldProducer
     process.PolyFit3DMagneticFieldProducer = PolyFit3DMagneticFieldProducer
     fieldlabel = "PolyFit3DMf"
     process.PolyFit3DMagneticFieldProducer.label = fieldlabel
+    process.geopro.MagneticFieldLabel = fieldlabel
     process.Geant4ePropagator.MagneticFieldLabel = fieldlabel
+    process.stripCPEESProducer.MagneticFieldLabel = fieldlabel
+    process.StripCPEfromTrackAngleESProducer.MagneticFieldLabel = fieldlabel
+    process.siPixelTemplateDBObjectESProducer.MagneticFieldLabel = fieldlabel
+    process.templates.MagneticFieldLabel = fieldlabel
+    process.trackrefit.MagneticFieldLabel = fieldlabel
 
     return process
 
diff --git a/RecoLocalTracker/SiPixelRecHits/interface/PixelCPETemplateRecoESProducer.h b/RecoLocalTracker/SiPixelRecHits/interface/PixelCPETemplateRecoESProducer.h
index 4a47426f44f23..2d0b8b73730de 100644
--- a/RecoLocalTracker/SiPixelRecHits/interface/PixelCPETemplateRecoESProducer.h
+++ b/RecoLocalTracker/SiPixelRecHits/interface/PixelCPETemplateRecoESProducer.h
@@ -15,6 +15,7 @@ class  PixelCPETemplateRecoESProducer: public edm::ESProducer{
  private:
   edm::ParameterSet pset_;
   bool DoLorentz_;
+  std::string fieldlabel_;
 };
 
 
diff --git a/RecoLocalTracker/SiPixelRecHits/plugins/PixelCPETemplateRecoESProducer.cc b/RecoLocalTracker/SiPixelRecHits/plugins/PixelCPETemplateRecoESProducer.cc
index 3bc2e3ff4bf57..f7d4a8a8ba824 100644
--- a/RecoLocalTracker/SiPixelRecHits/plugins/PixelCPETemplateRecoESProducer.cc
+++ b/RecoLocalTracker/SiPixelRecHits/plugins/PixelCPETemplateRecoESProducer.cc
@@ -27,6 +27,7 @@ PixelCPETemplateRecoESProducer::PixelCPETemplateRecoESProducer(const edm::Parame
   DoLorentz_ = p.existsAs<bool>("DoLorentz")?p.getParameter<bool>("DoLorentz"):false;
 
   pset_ = p;
+  fieldlabel_ = p.getParameter<std::string>("MagneticFieldLabel");
   setWhatProduced(this,myname);
 
   //std::cout<<" from ES Producer Templates "<<myname<<" "<<DoLorentz_<<std::endl;  //dk
@@ -39,7 +40,7 @@ std::unique_ptr<PixelClusterParameterEstimator>
 PixelCPETemplateRecoESProducer::produce(const TkPixelCPERecord & iRecord){ 
 
   ESHandle<MagneticField> magfield;
-  iRecord.getRecord<IdealMagneticFieldRecord>().get(magfield );
+  iRecord.getRecord<IdealMagneticFieldRecord>().get(fieldlabel_, magfield );
 
   edm::ESHandle<TrackerGeometry> pDD;
   iRecord.getRecord<TrackerDigiGeometryRecord>().get( pDD );
diff --git a/RecoLocalTracker/SiPixelRecHits/python/PixelCPETemplateReco_cfi.py b/RecoLocalTracker/SiPixelRecHits/python/PixelCPETemplateReco_cfi.py
index acebe3320de40..b358f448229da 100644
--- a/RecoLocalTracker/SiPixelRecHits/python/PixelCPETemplateReco_cfi.py
+++ b/RecoLocalTracker/SiPixelRecHits/python/PixelCPETemplateReco_cfi.py
@@ -15,7 +15,8 @@
     # True in Run II for offline RECO
     DoLorentz = cms.bool(True),
  
-    LoadTemplatesFromDB = cms.bool(True)
+    LoadTemplatesFromDB = cms.bool(True),
+    MagneticFieldLabel = cms.string(""),
 
 )
 
diff --git a/RecoLocalTracker/SiStripRecHitConverter/plugins/StripCPEESProducer.cc b/RecoLocalTracker/SiStripRecHitConverter/plugins/StripCPEESProducer.cc
index 0f5f82710034e..40c099e6074d1 100644
--- a/RecoLocalTracker/SiStripRecHitConverter/plugins/StripCPEESProducer.cc
+++ b/RecoLocalTracker/SiStripRecHitConverter/plugins/StripCPEESProducer.cc
@@ -30,6 +30,7 @@ StripCPEESProducer::StripCPEESProducer(const edm::ParameterSet & p)
 
   cpeNum = enumMap[type];
   parametersPSet = (p.exists("parameters") ? p.getParameter<edm::ParameterSet>("parameters") : p);
+  fieldlabel = p.getParameter<std::string>("MagneticFieldLabel");
   setWhatProduced(this,name);
 }
 
@@ -37,7 +38,7 @@ std::unique_ptr<StripClusterParameterEstimator> StripCPEESProducer::
 produce(const TkStripCPERecord & iRecord) 
 { 
   edm::ESHandle<TrackerGeometry> pDD;  iRecord.getRecord<TrackerDigiGeometryRecord>().get( pDD );
-  edm::ESHandle<MagneticField> magfield;  iRecord.getRecord<IdealMagneticFieldRecord>().get(magfield );
+  edm::ESHandle<MagneticField> magfield;  iRecord.getRecord<IdealMagneticFieldRecord>().get(fieldlabel, magfield );
   edm::ESHandle<SiStripLorentzAngle> lorentzAngle;   iRecord.getRecord<SiStripLorentzAngleDepRcd>().get(lorentzAngle);
   edm::ESHandle<SiStripBackPlaneCorrection> backPlaneCorrection; iRecord.getRecord<SiStripBackPlaneCorrectionDepRcd>().get(backPlaneCorrection);
   edm::ESHandle<SiStripConfObject> confObj;  iRecord.getRecord<SiStripConfObjectRcd>().get(confObj);
diff --git a/RecoLocalTracker/SiStripRecHitConverter/plugins/StripCPEESProducer.h b/RecoLocalTracker/SiStripRecHitConverter/plugins/StripCPEESProducer.h
index f7c0ff1de71b0..df4f86116d92b 100644
--- a/RecoLocalTracker/SiStripRecHitConverter/plugins/StripCPEESProducer.h
+++ b/RecoLocalTracker/SiStripRecHitConverter/plugins/StripCPEESProducer.h
@@ -25,6 +25,7 @@ class  StripCPEESProducer: public edm::ESProducer {
 
   CPE_t cpeNum;
   edm::ParameterSet parametersPSet;
+  std::string fieldlabel;
 
 };
 #endif
diff --git a/RecoLocalTracker/SiStripRecHitConverter/python/StripCPEESProducer_cfi.py b/RecoLocalTracker/SiStripRecHitConverter/python/StripCPEESProducer_cfi.py
index c3dc44c778d9d..3190d39f4efc1 100644
--- a/RecoLocalTracker/SiStripRecHitConverter/python/StripCPEESProducer_cfi.py
+++ b/RecoLocalTracker/SiStripRecHitConverter/python/StripCPEESProducer_cfi.py
@@ -3,5 +3,6 @@
 stripCPEESProducer = cms.ESProducer("StripCPEESProducer",
                                     ComponentName = cms.string('stripCPE'),
                                     ComponentType = cms.string('SimpleStripCPE'),
-                                    parameters    = cms.PSet()
+                                    parameters    = cms.PSet(),
+                                    MagneticFieldLabel = cms.string(""),
 )
diff --git a/SimG4Core/GeometryProducer/interface/GeometryProducer.h b/SimG4Core/GeometryProducer/interface/GeometryProducer.h
index df799048a03dc..3938c7d7be1f8 100644
--- a/SimG4Core/GeometryProducer/interface/GeometryProducer.h
+++ b/SimG4Core/GeometryProducer/interface/GeometryProducer.h
@@ -47,6 +47,7 @@ class GeometryProducer : public edm::one::EDProducer<edm::one::SharedResources,
   void updateMagneticField(edm::EventSetup const &es);
 
   G4RunManagerKernel *m_kernel;
+  std::string m_fieldlabel;
   edm::ParameterSet m_pField;
   SimActivityRegistry m_registry;
   std::vector<std::shared_ptr<SimWatcher>> m_watchers;
diff --git a/SimG4Core/GeometryProducer/src/GeometryProducer.cc b/SimG4Core/GeometryProducer/src/GeometryProducer.cc
index fbd6dc4d238cb..39d7a817c7d6d 100644
--- a/SimG4Core/GeometryProducer/src/GeometryProducer.cc
+++ b/SimG4Core/GeometryProducer/src/GeometryProducer.cc
@@ -59,6 +59,7 @@ static void createWatchers(const edm::ParameterSet &iP,
 
 GeometryProducer::GeometryProducer(edm::ParameterSet const &p)
     : m_kernel(nullptr),
+      m_fieldlabel(p.getParameter<std::string>("MagneticFieldLabel")),
       m_pField(p.getParameter<edm::ParameterSet>("MagneticField")),
       m_attach(nullptr),
       m_p(p),
@@ -84,7 +85,7 @@ void GeometryProducer::updateMagneticField(edm::EventSetup const &es) {
   if (m_pUseMagneticField) {
     // setup the magnetic field
     edm::ESHandle<MagneticField> pMF;
-    es.get<IdealMagneticFieldRecord>().get(pMF);
+    es.get<IdealMagneticFieldRecord>().get(m_fieldlabel, pMF);
     const GlobalPoint g(0., 0., 0.);
     edm::LogInfo("GeometryProducer") << "B-field(T) at (0,0,0)(cm): " << pMF->inTesla(g);
 
diff --git a/TrackPropagation/Geant4e/python/geantRefit_cff.py b/TrackPropagation/Geant4e/python/geantRefit_cff.py
index 2463f49a09564..4689ab19993b0 100644
--- a/TrackPropagation/Geant4e/python/geantRefit_cff.py
+++ b/TrackPropagation/Geant4e/python/geantRefit_cff.py
@@ -10,6 +10,7 @@
 geopro = cms.EDProducer("GeometryProducer",
      UseMagneticField = cms.bool(True),
      UseSensitiveDetectors = cms.bool(False),
+     MagneticFieldLabel = cms.string(""),
      MagneticField = cms.PSet(
          UseLocalMagFieldManager = cms.bool(False),
          Verbosity = cms.untracked.bool(False),
diff --git a/TrackingTools/TransientTrack/plugins/TransientTrackBuilderESProducer.cc b/TrackingTools/TransientTrack/plugins/TransientTrackBuilderESProducer.cc
index ff7518b73ce6e..9866186d151a1 100644
--- a/TrackingTools/TransientTrack/plugins/TransientTrackBuilderESProducer.cc
+++ b/TrackingTools/TransientTrack/plugins/TransientTrackBuilderESProducer.cc
@@ -15,6 +15,7 @@ TransientTrackBuilderESProducer::TransientTrackBuilderESProducer(const edm::Para
 {
   std::string myname = p.getParameter<std::string>("ComponentName");
   pset_ = p;
+  fieldlabel_ = p.getParameter<std::string>("MagneticFieldLabel");
   setWhatProduced(this,myname);
 }
 
@@ -24,7 +25,7 @@ std::unique_ptr<TransientTrackBuilder>
 TransientTrackBuilderESProducer::produce(const TransientTrackRecord & iRecord){ 
 
   edm::ESHandle<MagneticField> magfield;
-  iRecord.getRecord<IdealMagneticFieldRecord>().get( magfield );     
+  iRecord.getRecord<IdealMagneticFieldRecord>().get(fieldlabel_, magfield );
   edm::ESHandle<GlobalTrackingGeometry> theTrackingGeometry;
   iRecord.getRecord<GlobalTrackingGeometryRecord>().get(theTrackingGeometry); 
 
diff --git a/TrackingTools/TransientTrack/plugins/TransientTrackBuilderESProducer.h b/TrackingTools/TransientTrack/plugins/TransientTrackBuilderESProducer.h
index 9a94953e3e8d3..eb751673b7571 100644
--- a/TrackingTools/TransientTrack/plugins/TransientTrackBuilderESProducer.h
+++ b/TrackingTools/TransientTrack/plugins/TransientTrackBuilderESProducer.h
@@ -17,6 +17,7 @@ class  TransientTrackBuilderESProducer: public edm::ESProducer{
   std::unique_ptr<TransientTrackBuilder> produce(const TransientTrackRecord &);
  private:
   edm::ParameterSet pset_;
+  std::string fieldlabel_;
 };
 
 
diff --git a/TrackingTools/TransientTrack/python/TransientTrackBuilder_cfi.py b/TrackingTools/TransientTrack/python/TransientTrackBuilder_cfi.py
index 47ee3465f719d..c599a372d5ccc 100644
--- a/TrackingTools/TransientTrack/python/TransientTrackBuilder_cfi.py
+++ b/TrackingTools/TransientTrack/python/TransientTrackBuilder_cfi.py
@@ -2,7 +2,8 @@
 
 
 TransientTrackBuilderESProducer = cms.ESProducer("TransientTrackBuilderESProducer",
-    ComponentName = cms.string('TransientTrackBuilder')
+    ComponentName = cms.string('TransientTrackBuilder'),
+    MagneticFieldLabel = cms.string(""),
 )