File indexing completed on 2024-04-06 12:25:28
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
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
0055
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 }