/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/RecoMuon/MuonIdentification/src/MuonKinkFinder.cc

Go to the documentation of this file.

#include "RecoMuon/MuonIdentification/interface/MuonKinkFinder.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "TrackingTools/PatternTools/interface/Trajectory.h"

MuonKinkFinder::MuonKinkFinder(const edm::ParameterSet &iConfig) :
    diagonalOnly_(iConfig.getParameter<bool>("diagonalOnly")),
    usePosition_(iConfig.getParameter<bool>("usePosition")),
    refitter_(iConfig)
{ 
}

MuonKinkFinder::~MuonKinkFinder() { }

void MuonKinkFinder::init(const edm::EventSetup &iSetup) { 
    refitter_.setServices(iSetup);
}

bool MuonKinkFinder::fillTrkKink(reco::MuonQuality & quality, const reco::Track &track) const {
    std::vector<Trajectory> traj  = refitter_.transform(track);
    if (traj.size() != 1) {
        quality.trkKink = 999;
        quality.tkKink_position = math::XYZPoint(0,0,0);
        return false;
    }
    return fillTrkKink(quality, traj.front());
}

bool MuonKinkFinder::fillTrkKink(reco::MuonQuality & quality, const Trajectory &trajectory) const {
    const std::vector<TrajectoryMeasurement> &tms = trajectory.measurements();
    quality.trkKink = -1.0;
    quality.tkKink_position = math::XYZPoint(0,0,0);
    bool found = false; 
    for (int itm = 3, nm = tms.size() - 3;  itm < nm; ++itm) {
        TrajectoryStateOnSurface pre  = tms[itm].forwardPredictedState();
        TrajectoryStateOnSurface post = tms[itm].backwardPredictedState();
        if (!pre.isValid() || !post.isValid()) continue;
        found = true;
        double c2f = getChi2(pre,post);
        if (c2f > quality.trkKink) { 
            quality.trkKink = c2f; 
            GlobalPoint pos = (tms[itm].updatedState().isValid() ? tms[itm].updatedState() : pre).globalPosition();
            quality.tkKink_position = math::XYZPoint(pos.x(), pos.y(), pos.z());
        }
    }
    if (!found) quality.trkKink = 999;
    return found;
}

double MuonKinkFinder::getChi2(const TrajectoryStateOnSurface &start, const TrajectoryStateOnSurface &other) const {
    if (!start.hasError() && !other.hasError()) throw cms::Exception("LogicError") << "At least one of the two states must have errors to make chi2s.\n";
    AlgebraicSymMatrix55 cov;
    if (start.hasError()) cov += start.localError().matrix();
    if (other.hasError()) cov += other.localError().matrix();
    cropAndInvert(cov);
    AlgebraicVector5 diff(start.localParameters().mixedFormatVector() - other.localParameters().mixedFormatVector());
    return ROOT::Math::Similarity(diff, cov);
}

void MuonKinkFinder::cropAndInvert(AlgebraicSymMatrix55 &cov) const {
    if (usePosition_) {
        if (diagonalOnly_) {
            for (size_t i = 0; i < 5; ++i) { for (size_t j = i+1; j < 5; ++j) {
                cov(i,j) = 0;
            } }
        }
        cov.Invert();
    } else {
        // get 3x3 covariance
        AlgebraicSymMatrix33 momCov = cov.Sub<AlgebraicSymMatrix33>(0,0); // get 3x3 matrix
        if (diagonalOnly_) { momCov(0,1) = 0; momCov(0,2) = 0; momCov(1,2) = 0; }
        // invert        
        momCov.Invert();
        // place it
        cov.Place_at(momCov,0,0);
        // zero the rest
        for (size_t i = 3; i < 5; ++i) { for (size_t j = i; j < 5; ++j) { 
            cov(i,j) = 0;
        } }
    }
}