File indexing completed on 2021-06-24 02:11:19
0001 #ifndef TrackVertexArbitration_H
0002 #define TrackVertexArbitration_H
0003 #include <memory>
0004 #include <set>
0005 #include <unordered_map>
0006
0007 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0008 #include "TrackingTools/TransientTrack/interface/CandidatePtrTransientTrack.h"
0009 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0010 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0011 #include "TrackingTools/IPTools/interface/IPTools.h"
0012
0013 #include "FWCore/Framework/interface/EDProducer.h"
0014 #include "FWCore/Framework/interface/Event.h"
0015 #include "FWCore/Framework/interface/MakerMacros.h"
0016 #include "FWCore/Utilities/interface/InputTag.h"
0017 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0018
0019 #include "DataFormats/Common/interface/Handle.h"
0020 #include "DataFormats/TrackReco/interface/Track.h"
0021 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0022 #include "DataFormats/VertexReco/interface/Vertex.h"
0023 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0024 #include "RecoVertex/VertexTools/interface/VertexDistance3D.h"
0025
0026 #include "DataFormats/TrackReco/interface/Track.h"
0027 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0028 #include "DataFormats/VertexReco/interface/Vertex.h"
0029 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0030 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0031
0032 #include "RecoVertex/AdaptiveVertexFit/interface/AdaptiveVertexFitter.h"
0033 #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexUpdator.h"
0034 #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexTrackCompatibilityEstimator.h"
0035 #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexSmoother.h"
0036 #include "TrackingTools/GeomPropagators/interface/AnalyticalImpactPointExtrapolator.h"
0037 #include "RecoVertex/VertexPrimitives/interface/ConvertToFromReco.h"
0038 #include "DataFormats/Candidate/interface/CandidateFwd.h"
0039 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0040
0041 #include "RecoVertex/AdaptiveVertexFinder/interface/SVTimeHelpers.h"
0042
0043 #include "FWCore/Utilities/interface/isFinite.h"
0044
0045
0046
0047
0048 namespace svhelper {
0049 inline double cov33(const reco::Vertex &sv) { return sv.covariance(3, 3); }
0050 inline double cov33(const reco::VertexCompositePtrCandidate &sv) { return sv.vertexCovariance(3, 3); }
0051 }
0052
0053 template <class VTX>
0054 class TrackVertexArbitration {
0055 public:
0056 TrackVertexArbitration(const edm::ParameterSet ¶ms);
0057
0058 std::vector<VTX> trackVertexArbitrator(edm::Handle<reco::BeamSpot> &beamSpot,
0059 const reco::Vertex &pv,
0060 std::vector<reco::TransientTrack> &selectedTracks,
0061 std::vector<VTX> &secondaryVertices);
0062
0063 private:
0064 bool trackFilterArbitrator(const reco::TransientTrack &track) const;
0065
0066 edm::InputTag primaryVertexCollection;
0067 edm::InputTag secondaryVertexCollection;
0068 edm::InputTag trackCollection;
0069 edm::InputTag beamSpotCollection;
0070 double dRCut;
0071 double distCut;
0072 double sigCut;
0073 double dLenFraction;
0074 double fitterSigmacut;
0075 double fitterTini;
0076 double fitterRatio;
0077 int trackMinLayers;
0078 double trackMinPt;
0079 int trackMinPixels;
0080 double maxTimeSignificance;
0081 double sign;
0082 };
0083
0084 #include "DataFormats/GeometryVector/interface/VectorUtil.h"
0085 template <class VTX>
0086 TrackVertexArbitration<VTX>::TrackVertexArbitration(const edm::ParameterSet ¶ms)
0087 : primaryVertexCollection(params.getParameter<edm::InputTag>("primaryVertices")),
0088 secondaryVertexCollection(params.getParameter<edm::InputTag>("secondaryVertices")),
0089 trackCollection(params.getParameter<edm::InputTag>("tracks")),
0090 beamSpotCollection(params.getParameter<edm::InputTag>("beamSpot")),
0091 dRCut(params.getParameter<double>("dRCut")),
0092 distCut(params.getParameter<double>("distCut")),
0093 sigCut(params.getParameter<double>("sigCut")),
0094 dLenFraction(params.getParameter<double>("dLenFraction")),
0095 fitterSigmacut(params.getParameter<double>("fitterSigmacut")),
0096 fitterTini(params.getParameter<double>("fitterTini")),
0097 fitterRatio(params.getParameter<double>("fitterRatio")),
0098 trackMinLayers(params.getParameter<int32_t>("trackMinLayers")),
0099 trackMinPt(params.getParameter<double>("trackMinPt")),
0100 trackMinPixels(params.getParameter<int32_t>("trackMinPixels")),
0101 maxTimeSignificance(params.getParameter<double>("maxTimeSignificance")) {
0102 sign = 1.0;
0103 if (dRCut < 0) {
0104 sign = -1.0;
0105 }
0106 dRCut *= dRCut;
0107 }
0108 template <class VTX>
0109 bool TrackVertexArbitration<VTX>::trackFilterArbitrator(const reco::TransientTrack &track) const {
0110 if (!track.isValid())
0111 return false;
0112 if (track.track().hitPattern().trackerLayersWithMeasurement() < trackMinLayers)
0113 return false;
0114 if (track.track().pt() < trackMinPt)
0115 return false;
0116 if (track.track().hitPattern().numberOfValidPixelHits() < trackMinPixels)
0117 return false;
0118
0119 return true;
0120 }
0121
0122 inline float trackWeight(const reco::Vertex &sv, const reco::TransientTrack &track) {
0123 return sv.trackWeight(track.trackBaseRef());
0124 }
0125 inline float trackWeight(const reco::VertexCompositePtrCandidate &sv, const reco::TransientTrack &tt) {
0126 const reco::CandidatePtrTransientTrack *cptt =
0127 dynamic_cast<const reco::CandidatePtrTransientTrack *>(tt.basicTransientTrack());
0128 if (cptt == nullptr)
0129 edm::LogError("DynamicCastingFailed") << "Casting of TransientTrack to CandidatePtrTransientTrack failed!";
0130 else {
0131 const reco::CandidatePtr &c = cptt->candidate();
0132 if (std::find(sv.daughterPtrVector().begin(), sv.daughterPtrVector().end(), c) != sv.daughterPtrVector().end())
0133 return 1.0;
0134 else
0135 return 0;
0136 }
0137 return 0;
0138 }
0139
0140 template <class VTX>
0141 std::vector<VTX> TrackVertexArbitration<VTX>::trackVertexArbitrator(edm::Handle<reco::BeamSpot> &beamSpot,
0142 const reco::Vertex &pv,
0143 std::vector<reco::TransientTrack> &selectedTracks,
0144 std::vector<VTX> &secondaryVertices) {
0145 using namespace reco;
0146
0147
0148 VertexDistance3D dist;
0149 AdaptiveVertexFitter theAdaptiveFitter(GeometricAnnealing(fitterSigmacut, fitterTini, fitterRatio),
0150 DefaultLinearizationPointFinder(),
0151 KalmanVertexUpdator<5>(),
0152 KalmanVertexTrackCompatibilityEstimator<5>(),
0153 KalmanVertexSmoother());
0154
0155 std::vector<VTX> recoVertices;
0156 VertexDistance3D vdist;
0157 GlobalPoint ppv(pv.position().x(), pv.position().y(), pv.position().z());
0158
0159 std::unordered_map<unsigned int, Measurement1D> cachedIP;
0160 for (typename std::vector<VTX>::const_iterator sv = secondaryVertices.begin(); sv != secondaryVertices.end(); ++sv) {
0161 const double svTime(sv->t()), svTimeCov(svhelper::cov33(*sv));
0162 const bool svHasTime = svTimeCov > 0.;
0163
0164 GlobalPoint ssv(sv->position().x(), sv->position().y(), sv->position().z());
0165 GlobalVector flightDir = ssv - ppv;
0166
0167 Measurement1D dlen =
0168 vdist.distance(pv, VertexState(RecoVertex::convertPos(sv->position()), RecoVertex::convertError(sv->error())));
0169 std::vector<reco::TransientTrack> selTracks;
0170 for (unsigned int itrack = 0; itrack < selectedTracks.size(); itrack++) {
0171 TransientTrack &tt = (selectedTracks)[itrack];
0172 if (!trackFilterArbitrator(tt))
0173 continue;
0174 tt.setBeamSpot(*beamSpot);
0175 float w = trackWeight(*sv, tt);
0176 Measurement1D ipv;
0177 if (cachedIP.count(itrack)) {
0178 ipv = cachedIP[itrack];
0179 } else {
0180 std::pair<bool, Measurement1D> ipvp = IPTools::absoluteImpactParameter3D(tt, pv);
0181 cachedIP[itrack] = ipvp.second;
0182 ipv = ipvp.second;
0183 }
0184
0185 AnalyticalImpactPointExtrapolator extrapolator(tt.field());
0186 TrajectoryStateOnSurface tsos =
0187 extrapolator.extrapolate(tt.impactPointState(), RecoVertex::convertPos(sv->position()));
0188 if (!tsos.isValid())
0189 continue;
0190 GlobalPoint refPoint = tsos.globalPosition();
0191 GlobalError refPointErr = tsos.cartesianError().position();
0192 Measurement1D isv =
0193 dist.distance(VertexState(RecoVertex::convertPos(sv->position()), RecoVertex::convertError(sv->error())),
0194 VertexState(refPoint, refPointErr));
0195
0196 float dR = Geom::deltaR2(((sign > 0) ? flightDir : -flightDir),
0197 tt.track());
0198
0199 double timeSig = 0.;
0200 if (svHasTime && edm::isFinite(tt.timeExt())) {
0201 double tError = std::sqrt(std::pow(tt.dtErrorExt(), 2) + svTimeCov);
0202 timeSig = std::abs(tt.timeExt() - svTime) / tError;
0203 }
0204
0205 if (w > 0 || (isv.significance() < sigCut && isv.value() < distCut && isv.value() < dlen.value() * dLenFraction &&
0206 timeSig < maxTimeSignificance)) {
0207 if ((isv.value() < ipv.value()) && isv.value() < distCut && isv.value() < dlen.value() * dLenFraction &&
0208 dR < dRCut && timeSig < maxTimeSignificance) {
0209 #ifdef VTXDEBUG
0210 if (w > 0.5)
0211 std::cout << " = ";
0212 else
0213 std::cout << " + ";
0214 #endif
0215 selTracks.push_back(tt);
0216 } else {
0217 #ifdef VTXDEBUG
0218 if (w > 0.5 && isv.value() > ipv.value())
0219 std::cout << " - ";
0220 else
0221 std::cout << " ";
0222 #endif
0223
0224 if (w > 0.5 && isv.value() <= ipv.value() && dR < dRCut && timeSig < maxTimeSignificance) {
0225 selTracks.push_back(tt);
0226 #ifdef VTXDEBUG
0227 std::cout << " = ";
0228 #endif
0229 }
0230 if (w > 0.5 && isv.value() <= ipv.value() && dR >= dRCut) {
0231 #ifdef VTXDEBUG
0232 std::cout << " - ";
0233 #endif
0234 }
0235 }
0236 #ifdef VTXDEBUG
0237
0238 std::cout << "t : " << itrack << " ref " << ref.key() << " w: " << w << " svip: " << isv.significance() << " "
0239 << isv.value() << " pvip: " << ipv.significance() << " " << ipv.value() << " res "
0240 << tt.track().residualX(0) << "," << tt.track().residualY(0)
0241
0242 #endif
0243
0244 } else {
0245 #ifdef VTXDEBUG
0246
0247 std::cout << " . t : " << itrack << " ref " << ref.key() << " w: " << w
0248 << " svip: " << isv.second.significance() << " " << isv.second.value()
0249 << " pvip: " << ipv.significance() << " " << ipv.value() << " dr: " << dR << std::endl;
0250 #endif
0251 }
0252 }
0253
0254 if (selTracks.size() >= 2) {
0255 TransientVertex singleFitVertex;
0256 singleFitVertex = theAdaptiveFitter.vertex(selTracks, ssv);
0257
0258 if (singleFitVertex.isValid()) {
0259 if (pv.covariance(3, 3) > 0.)
0260 svhelper::updateVertexTime(singleFitVertex);
0261 recoVertices.push_back(VTX(singleFitVertex));
0262
0263 #ifdef VTXDEBUG
0264 const VTX &extVertex = recoVertices.back();
0265 GlobalVector vtxDir = GlobalVector(extVertex.p4().X(), extVertex.p4().Y(), extVertex.p4().Z());
0266
0267 std::cout << " pointing : " << Geom::deltaR(extVertex.position() - pv.position(), vtxDir) << std::endl;
0268 #endif
0269 }
0270 }
0271 }
0272 return recoVertices;
0273 }
0274
0275 #endif