Merge pull request #28753 from cericeci/fixPatMuondXY

Fixed sign in pat::Muon::dB() to match one obtained from the best track definition
cms-sw · Feb 7, 2020 · 736475d · 736475d
2 parents 1edbed3 + 8dd5d67
commit 736475d
Show file tree

Hide file tree

Showing 3 changed files with 34 additions and 13 deletions.
diff --git a/DataFormats/TrackReco/interface/TrackBase.h b/DataFormats/TrackReco/interface/TrackBase.h
@@ -335,6 +335,12 @@ namespace reco {
     /// error on beta
     double betaError() const;
 
+    /// error on dxy with respect to a user-given reference point + uncertainty (i.e. reco::Vertex position)
+    double dxyError(Point const &vtx, math::Error<3>::type const &vertexCov) const;
+
+    /// error on dxy with respect to a user-given beamspot
+    double dxyError(const BeamSpot &theBeamSpot) const;
+
     /// fill SMatrix
     CovarianceMatrix &fill(CovarianceMatrix &v) const;
 
@@ -736,6 +742,11 @@ namespace reco {
   // error on beta
   inline double TrackBase::betaError() const { return std::sqrt(covbetabeta_); }
 
+  // error on dxy with respect to a given beamspot
+  inline double TrackBase::dxyError(const BeamSpot &theBeamSpot) const {
+    return dxyError(theBeamSpot.position(vz()), theBeamSpot.rotatedCovariance3D());
+  }
+
   // number of valid hits found
   inline unsigned short TrackBase::numberOfValidHits() const { return hitPattern_.numberOfValidHits(); }
 

diff --git a/DataFormats/TrackReco/src/TrackBase.cc b/DataFormats/TrackReco/src/TrackBase.cc
@@ -144,3 +144,16 @@ TrackBase::TrackAlgorithm TrackBase::algoByName(const std::string &name) {
   // cast
   return TrackAlgorithm(index);
 }
+
+double TrackBase::dxyError(Point const &vtx, math::Error<3>::type const &vertexCov) const {
+  // Gradient of TrackBase::dxy(const Point &myBeamSpot) with respect to track parameters. Using unrolled expressions to avoid calling for higher dimension matrices
+  // ( 0, 0, x_vert * cos(phi) + y_vert * sin(phi), 1, 0 )
+  // Gradient with respect to point parameters
+  // ( sin(phi), -cos(phi))
+  // Propagate covariance assuming cross-terms of the covariance between track and vertex parameters are 0
+  return std::sqrt((vtx.x() * px() + vtx.y() * py()) * (vtx.x() * px() + vtx.y() * py()) / (pt() * pt()) *
+                       covariance(i_phi, i_phi) +
+                   2 * (vtx.x() * px() + vtx.y() * py()) / pt() * covariance(i_phi, i_dxy) + covariance(i_dxy, i_dxy) +
+                   py() * py() / (pt() * pt()) * vertexCov(0, 0) - 2 * py() * px() / (pt() * pt()) * vertexCov(0, 1) +
+                   px() * px() / (pt() * pt()) * vertexCov(1, 1));
+}
diff --git a/PhysicsTools/PatAlgos/plugins/PATMuonProducer.cc b/PhysicsTools/PatAlgos/plugins/PATMuonProducer.cc
@@ -1097,27 +1097,24 @@ void PATMuonProducer::embedHighLevel(pat::Muon& aMuon,
   // Correct to PV
 
   // PV2D
+  aMuon.setDB(track->dxy(primaryVertex.position()),
+              track->dxyError(primaryVertex.position(), primaryVertex.covariance()),
+              pat::Muon::PV2D);
+
+  // PV3D
   std::pair<bool, Measurement1D> result =
-      IPTools::signedTransverseImpactParameter(tt, GlobalVector(track->px(), track->py(), track->pz()), primaryVertex);
+      IPTools::signedImpactParameter3D(tt, GlobalVector(track->px(), track->py(), track->pz()), primaryVertex);
   double d0_corr = result.second.value();
   double d0_err = primaryVertexIsValid ? result.second.error() : -1.0;
-  aMuon.setDB(d0_corr, d0_err, pat::Muon::PV2D);
-
-  // PV3D
-  result = IPTools::signedImpactParameter3D(tt, GlobalVector(track->px(), track->py(), track->pz()), primaryVertex);
-  d0_corr = result.second.value();
-  d0_err = primaryVertexIsValid ? result.second.error() : -1.0;
   aMuon.setDB(d0_corr, d0_err, pat::Muon::PV3D);
 
   // Correct to beam spot
-  // make a fake vertex out of beam spot
-  reco::Vertex vBeamspot(beamspot.position(), beamspot.rotatedCovariance3D());
 
   // BS2D
-  result = IPTools::signedTransverseImpactParameter(tt, GlobalVector(track->px(), track->py(), track->pz()), vBeamspot);
-  d0_corr = result.second.value();
-  d0_err = beamspotIsValid ? result.second.error() : -1.0;
-  aMuon.setDB(d0_corr, d0_err, pat::Muon::BS2D);
+  aMuon.setDB(track->dxy(beamspot), track->dxyError(beamspot), pat::Muon::BS2D);
+
+  // make a fake vertex out of beam spot
+  reco::Vertex vBeamspot(beamspot.position(), beamspot.rotatedCovariance3D());
 
   // BS3D
   result = IPTools::signedImpactParameter3D(tt, GlobalVector(track->px(), track->py(), track->pz()), vBeamspot);