MuonIdentification/src/MuonKinkFinder.cc

0001 #include "RecoMuon/MuonIdentification/interface/MuonKinkFinder.h"
0002 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0003 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0004
0005 MuonKinkFinder::MuonKinkFinder(const edm::ParameterSet &iConfig, edm::ConsumesCollector &iC)
0006     : diagonalOnly_(iConfig.getParameter<bool>("diagonalOnly")),
0007       usePosition_(iConfig.getParameter<bool>("usePosition")),
0008       refitter_(iConfig, iC) {}
0009
0010 MuonKinkFinder::~MuonKinkFinder() {}
0011
0012 void MuonKinkFinder::init(const edm::EventSetup &iSetup) { refitter_.setServices(iSetup); }
0013
0014 bool MuonKinkFinder::fillTrkKink(reco::MuonQuality &quality, const reco::Track &track) const {
0015   std::vector<Trajectory> traj = refitter_.transform(track);
0016   if (traj.size() != 1) {
0017     quality.trkKink = 999;
0018     quality.tkKink_position = math::XYZPoint(0, 0, 0);
0019     return false;
0020   }
0021   return fillTrkKink(quality, traj.front());
0022 }
0023
0024 bool MuonKinkFinder::fillTrkKink(reco::MuonQuality &quality, const Trajectory &trajectory) const {
0025   const std::vector<TrajectoryMeasurement> &tms = trajectory.measurements();
0026   quality.trkKink = -1.0;
0027   quality.tkKink_position = math::XYZPoint(0, 0, 0);
0028   bool found = false;
0029   for (int itm = 3, nm = tms.size() - 3; itm < nm; ++itm) {
0030     TrajectoryStateOnSurface pre = tms[itm].forwardPredictedState();
0031     TrajectoryStateOnSurface post = tms[itm].backwardPredictedState();
0032     if (!pre.isValid() || !post.isValid())
0033       continue;
0034     found = true;
0035     double c2f = getChi2(pre, post);
0036     if (c2f > quality.trkKink) {
0037       quality.trkKink = c2f;
0038       GlobalPoint pos = (tms[itm].updatedState().isValid() ? tms[itm].updatedState() : pre).globalPosition();
0039       quality.tkKink_position = math::XYZPoint(pos.x(), pos.y(), pos.z());
0040     }
0041   }
0042   if (!found)
0043     quality.trkKink = 999;
0044   return found;
0045 }
0046
0047 double MuonKinkFinder::getChi2(const TrajectoryStateOnSurface &start, const TrajectoryStateOnSurface &other) const {
0048   if (!start.hasError() && !other.hasError())
0049     throw cms::Exception("LogicError") << "At least one of the two states must have errors to make chi2s.\n";
0050   AlgebraicSymMatrix55 cov;
0051   if (start.hasError())
0052     cov += start.localError().matrix();
0053   if (other.hasError())
0054     cov += other.localError().matrix();
0055   cropAndInvert(cov);
0056   AlgebraicVector5 diff(start.localParameters().mixedFormatVector() - other.localParameters().mixedFormatVector());
0057   return ROOT::Math::Similarity(diff, cov);
0058 }
0059
0060 void MuonKinkFinder::cropAndInvert(AlgebraicSymMatrix55 &cov) const {
0061   if (usePosition_) {
0062     if (diagonalOnly_) {
0063       for (size_t i = 0; i < 5; ++i) {
0064         for (size_t j = i + 1; j < 5; ++j) {
0065           cov(i, j) = 0;
0066         }
0067       }
0068     }
0069     cov.Invert();
0070   } else {
0071     // get 3x3 covariance
0072     AlgebraicSymMatrix33 momCov = cov.Sub<AlgebraicSymMatrix33>(0, 0);  // get 3x3 matrix
0073     if (diagonalOnly_) {
0074       momCov(0, 1) = 0;
0075       momCov(0, 2) = 0;
0076       momCov(1, 2) = 0;
0077     }
0078     // invert
0079     momCov.Invert();
0080     // place it
0081     cov.Place_at(momCov, 0, 0);
0082     // zero the rest
0083     for (size_t i = 3; i < 5; ++i) {
0084       for (size_t j = i; j < 5; ++j) {
0085         cov(i, j) = 0;
0086       }
0087     }
0088   }
0089 }