JetSignalVertexCompatibilityAlgo.cc

Go to the documentation of this file.

#include <functional>
#include <cmath>
#include <map>

#include "DataFormats/Common/interface/RefToBase.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/VertexReco/interface/Vertex.h"
#include "DataFormats/VertexReco/interface/VertexFwd.h"

#include "TrackingTools/GeomPropagators/interface/AnalyticalImpactPointExtrapolator.h"
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
#include "TrackingTools/TransientTrack/interface/TransientTrack.h"
#include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
#include "TrackingTools/Records/interface/TransientTrackRecord.h"

#include "RecoVertex/VertexPrimitives/interface/ConvertToFromReco.h"

#include "RecoJets/JetAssociationAlgorithms/interface/JetSignalVertexCompatibilityAlgo.h"

using namespace reco;

// helper
template<typename T>
inline bool
JetSignalVertexCompatibilityAlgo::RefToBaseLess<T>::operator()(
    const edm::RefToBase<T> &r1, const edm::RefToBase<T> &r2) const
{ return r1.id() < r2.id() || (r1.id() == r2.id() && r1.key() < r2.key()); }


JetSignalVertexCompatibilityAlgo::JetSignalVertexCompatibilityAlgo(
                    double cut, double temperature) :
    cut(cut), temperature(temperature)
{}

JetSignalVertexCompatibilityAlgo::~JetSignalVertexCompatibilityAlgo()
{}

double JetSignalVertexCompatibilityAlgo::trackVertexCompat(
            const Vertex &vtx, const TransientTrack &track)
{
    GlobalPoint point1 = RecoVertex::convertPos(vtx.position());
    AnalyticalImpactPointExtrapolator extrap(track.field());
    TrajectoryStateOnSurface tsos =
            extrap.extrapolate(track.impactPointState(), point1);

    if (!tsos.isValid())
        return 1.0e6;

    GlobalPoint point2 = tsos.globalPosition();
    ROOT::Math::SVector<double, 3> dir(point1.x() - point2.x(),
                                       point1.y() - point2.y(),
                                       point1.z() - point2.z());
    GlobalError cov = RecoVertex::convertError(vtx.covariance()) +
                      tsos.cartesianError().position();

    return ROOT::Math::Mag2(dir) /
           std::sqrt(ROOT::Math::Similarity(cov.matrix(), dir));
}

const TransientTrack&
JetSignalVertexCompatibilityAlgo::convert(const TrackBaseRef &track) const
{
    TransientTrackMap::iterator pos = trackMap.lower_bound(track);
    if (pos != trackMap.end() && pos->first == track)
        return pos->second;

    // the castTo will only work with regular, i.e. no GsfTracks
    // the interface is not intrinsically polymorph...
    return trackMap.insert(pos,
        std::make_pair(track, trackBuilder->build(
                    track.castTo<TrackRef>())))->second;
}

double JetSignalVertexCompatibilityAlgo::activation(double compat) const
{
    return 1. / (std::exp((compat - cut) / temperature) + 1.);
}

std::vector<float>
JetSignalVertexCompatibilityAlgo::compatibility(
                const reco::VertexCollection &vertices,
                const reco::TrackRefVector &tracks) const
{
    std::vector<float> result(vertices.size(), 0.);
    float sum = 0.;

    for(TrackRefVector::const_iterator track = tracks.begin();
        track != tracks.end(); ++track) {
        const TransientTrack &transientTrack =
                    convert(TrackBaseRef(*track));

        for(unsigned int i = 0; i < vertices.size(); i++) {
            double compat =
                trackVertexCompat(vertices[i], transientTrack);
            double contribution =
                activation(compat) * (*track)->pt();

            result[i] += contribution;
            sum += contribution;
        }
    }

    if (sum < 1.0e-9) {
        for(unsigned int i = 0; i < result.size(); i++)
            result[i] = 1.0 / result.size();
    } else {
        for(unsigned int i = 0; i < result.size(); i++)
            result[i] /= sum;
    }

    return result;
}

void JetSignalVertexCompatibilityAlgo::resetEvent(
                    const TransientTrackBuilder *builder)
{
    trackMap.clear();
    trackBuilder = builder;
}