JetAssociationAlgorithms/src/JetSignalVertexCompatibilityAlgo.cc

0001 #include <functional>
0002 #include <cmath>
0003 #include <map>
0004
0005 #include "DataFormats/Common/interface/RefToBase.h"
0006 #include "DataFormats/TrackReco/interface/Track.h"
0007 #include "DataFormats/VertexReco/interface/Vertex.h"
0008 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0009
0010 #include "TrackingTools/GeomPropagators/interface/AnalyticalImpactPointExtrapolator.h"
0011 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0012 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0013 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0014 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0015
0016 #include "RecoVertex/VertexPrimitives/interface/ConvertToFromReco.h"
0017
0018 #include "RecoJets/JetAssociationAlgorithms/interface/JetSignalVertexCompatibilityAlgo.h"
0019
0020 using namespace reco;
0021
0022 // helper
0023 template <typename T>
0024 inline bool JetSignalVertexCompatibilityAlgo::RefToBaseLess<T>::operator()(const edm::RefToBase<T> &r1,
0025                                                                            const edm::RefToBase<T> &r2) const {
0026   return r1.id() < r2.id() || (r1.id() == r2.id() && r1.key() < r2.key());
0027 }
0028
0029 JetSignalVertexCompatibilityAlgo::JetSignalVertexCompatibilityAlgo(double cut, double temperature)
0030     : cut(cut), temperature(temperature) {}
0031
0032 JetSignalVertexCompatibilityAlgo::~JetSignalVertexCompatibilityAlgo() {}
0033
0034 double JetSignalVertexCompatibilityAlgo::trackVertexCompat(const Vertex &vtx, const TransientTrack &track) {
0035   GlobalPoint point1 = RecoVertex::convertPos(vtx.position());
0036   AnalyticalImpactPointExtrapolator extrap(track.field());
0037   TrajectoryStateOnSurface tsos = extrap.extrapolate(track.impactPointState(), point1);
0038
0039   if (!tsos.isValid())
0040     return 1.0e6;
0041
0042   GlobalPoint point2 = tsos.globalPosition();
0043   ROOT::Math::SVector<double, 3> dir(point1.x() - point2.x(), point1.y() - point2.y(), point1.z() - point2.z());
0044   GlobalError cov = RecoVertex::convertError(vtx.covariance()) + tsos.cartesianError().position();
0045
0046   return ROOT::Math::Mag2(dir) / std::sqrt(ROOT::Math::Similarity(cov.matrix(), dir));
0047 }
0048
0049 const TransientTrack &JetSignalVertexCompatibilityAlgo::convert(const TrackBaseRef &track) const {
0050   TransientTrackMap::iterator pos = trackMap.lower_bound(track);
0051   if (pos != trackMap.end() && pos->first == track)
0052     return pos->second;
0053
0054   // the castTo will only work with regular, i.e. no GsfTracks
0055   // the interface is not intrinsically polymorph...
0056   return trackMap.insert(pos, std::make_pair(track, trackBuilder->build(track.castTo<TrackRef>())))->second;
0057 }
0058
0059 double JetSignalVertexCompatibilityAlgo::activation(double compat) const {
0060   return 1. / (std::exp((compat - cut) / temperature) + 1.);
0061 }
0062
0063 std::vector<float> JetSignalVertexCompatibilityAlgo::compatibility(const reco::VertexCollection &vertices,
0064                                                                    const reco::TrackRefVector &tracks) const {
0065   std::vector<float> result(vertices.size(), 0.);
0066   float sum = 0.;
0067
0068   for (TrackRefVector::const_iterator track = tracks.begin(); track != tracks.end(); ++track) {
0069     const TransientTrack &transientTrack = convert(TrackBaseRef(*track));
0070
0071     for (unsigned int i = 0; i < vertices.size(); i++) {
0072       double compat = trackVertexCompat(vertices[i], transientTrack);
0073       double contribution = activation(compat) * (*track)->pt();
0074
0075       result[i] += contribution;
0076       sum += contribution;
0077     }
0078   }
0079
0080   if (sum < 1.0e-9) {
0081     for (unsigned int i = 0; i < result.size(); i++)
0082       result[i] = 1.0 / result.size();
0083   } else {
0084     for (unsigned int i = 0; i < result.size(); i++)
0085       result[i] /= sum;
0086   }
0087
0088   return result;
0089 }
0090
0091 void JetSignalVertexCompatibilityAlgo::resetEvent(const TransientTrackBuilder *builder) {
0092   trackMap.clear();
0093   trackBuilder = builder;
0094 }