Merge pull request cms-sw#1 from WMass/WmassNanoProd_10_6_26

Muon calibration updates

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;
 

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_;
 

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,13 +312,17 @@ 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");
 
     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,27 +2221,44 @@ 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;
               }
               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;

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

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();

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(""),
+                                                   )
 

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"/>

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();
   }
 }

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;
 }
 

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);
 }
 

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",

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