File indexing completed on 2024-04-06 12:29:02
0001 #include "RecoVertex/AdaptiveVertexFinder/interface/AdaptiveVertexReconstructor.h"
0002 #include "RecoVertex/VertexTools/interface/GeometricAnnealing.h"
0003 #include "FWCore/Utilities/interface/EDMException.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005 #include "RecoVertex/VertexTools/interface/DummyVertexSmoother.h"
0006 #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexSmoother.h"
0007 #include <algorithm>
0008
0009 using namespace std;
0010
0011 TransientVertex AdaptiveVertexReconstructor::cleanUp(const TransientVertex& old) const {
0012 vector<reco::TransientTrack> trks = old.originalTracks();
0013 vector<reco::TransientTrack> newtrks;
0014 TransientVertex::TransientTrackToFloatMap mp;
0015 static const float minweight = 1.e-8;
0016 for (vector<reco::TransientTrack>::const_iterator i = trks.begin(); i != trks.end(); ++i) {
0017 if (old.trackWeight(*i) > minweight) {
0018 newtrks.push_back(*i);
0019 mp[*i] = old.trackWeight(*i);
0020 }
0021 }
0022
0023 TransientVertex ret;
0024
0025 if (old.hasPrior()) {
0026 VertexState priorstate(old.priorPosition(), old.priorError());
0027 ret = TransientVertex(priorstate, old.vertexState(), newtrks, old.totalChiSquared(), old.degreesOfFreedom());
0028 } else {
0029 ret = TransientVertex(old.vertexState(), newtrks, old.totalChiSquared(), old.degreesOfFreedom());
0030 }
0031 ret.weightMap(mp);
0032
0033 if (old.hasRefittedTracks()) {
0034
0035 vector<reco::TransientTrack> newrfs;
0036 vector<reco::TransientTrack> oldrfs = old.refittedTracks();
0037 vector<reco::TransientTrack>::const_iterator origtrkiter = trks.begin();
0038 for (vector<reco::TransientTrack>::const_iterator i = oldrfs.begin(); i != oldrfs.end(); ++i) {
0039 if (old.trackWeight(*origtrkiter) > minweight) {
0040 newrfs.push_back(*i);
0041 }
0042 origtrkiter++;
0043 }
0044 if (!newrfs.empty())
0045 ret.refittedTracks(newrfs);
0046 }
0047
0048 if (ret.refittedTracks().size() > ret.originalTracks().size()) {
0049 edm::LogError("AdaptiveVertexReconstructor") << "More refitted tracks than original tracks!";
0050 }
0051
0052 return ret;
0053 }
0054
0055 void AdaptiveVertexReconstructor::erase(const TransientVertex& newvtx,
0056 set<reco::TransientTrack>& remainingtrks,
0057 float w) const {
0058
0059
0060
0061
0062 const vector<reco::TransientTrack>& origtrks = newvtx.originalTracks();
0063 for (vector<reco::TransientTrack>::const_iterator i = origtrks.begin(); i != origtrks.end(); ++i) {
0064 double weight = newvtx.trackWeight(*i);
0065 if (weight > w) {
0066 remainingtrks.erase(*i);
0067 };
0068 };
0069 }
0070
0071 AdaptiveVertexReconstructor::AdaptiveVertexReconstructor(float primcut, float seccut, float min_weight, bool smoothing)
0072 : thePrimaryFitter(nullptr), theSecondaryFitter(nullptr), theMinWeight(min_weight), theWeightThreshold(0.001) {
0073 setupFitters(primcut, 256., 0.25, seccut, 256., 0.25, smoothing);
0074 }
0075
0076 AdaptiveVertexReconstructor::~AdaptiveVertexReconstructor() {
0077 if (thePrimaryFitter)
0078 delete thePrimaryFitter;
0079 if (theSecondaryFitter)
0080 delete theSecondaryFitter;
0081 }
0082
0083 void AdaptiveVertexReconstructor::setupFitters(
0084 float primcut, float primT, float primr, float seccut, float secT, float secr, bool smoothing) {
0085 VertexSmoother<5>* smoother;
0086 if (smoothing) {
0087 smoother = new KalmanVertexSmoother();
0088 } else {
0089 smoother = new DummyVertexSmoother<5>();
0090 }
0091
0092 if (thePrimaryFitter)
0093 delete thePrimaryFitter;
0094 if (theSecondaryFitter)
0095 delete theSecondaryFitter;
0096
0097
0098
0099
0100
0101 thePrimaryFitter = new AdaptiveVertexFitter(GeometricAnnealing(primcut, primT, primr),
0102 DefaultLinearizationPointFinder(),
0103 KalmanVertexUpdator<5>(),
0104 KalmanVertexTrackCompatibilityEstimator<5>(),
0105 *smoother);
0106 thePrimaryFitter->setWeightThreshold(theWeightThreshold);
0107
0108 theSecondaryFitter = new AdaptiveVertexFitter(GeometricAnnealing(seccut, secT, secr),
0109 DefaultLinearizationPointFinder(),
0110 KalmanVertexUpdator<5>(),
0111 KalmanVertexTrackCompatibilityEstimator<5>(),
0112 *smoother);
0113 theSecondaryFitter->setWeightThreshold(0.);
0114
0115
0116
0117 delete smoother;
0118 }
0119
0120 AdaptiveVertexReconstructor::AdaptiveVertexReconstructor(const edm::ParameterSet& m)
0121 : thePrimaryFitter(nullptr), theSecondaryFitter(nullptr), theMinWeight(0.5), theWeightThreshold(0.001) {
0122 float primcut = 2.0;
0123 float seccut = 6.0;
0124 bool smoothing = false;
0125
0126
0127 float primT = 256.;
0128 float primr = 0.25;
0129 float secT = 256.;
0130 float secr = 0.25;
0131
0132 try {
0133 primcut = m.getParameter<double>("primcut");
0134 primT = m.getParameter<double>("primT");
0135 primr = m.getParameter<double>("primr");
0136 seccut = m.getParameter<double>("seccut");
0137 secT = m.getParameter<double>("secT");
0138 secr = m.getParameter<double>("secr");
0139 theMinWeight = m.getParameter<double>("minweight");
0140 theWeightThreshold = m.getParameter<double>("weightthreshold");
0141 smoothing = m.getParameter<bool>("smoothing");
0142 } catch (edm::Exception& e) {
0143 edm::LogError("AdaptiveVertexReconstructor") << e.what();
0144 }
0145
0146 setupFitters(primcut, primT, primr, seccut, secT, secr, smoothing);
0147 }
0148
0149 vector<TransientVertex> AdaptiveVertexReconstructor::vertices(const vector<reco::TransientTrack>& t,
0150 const reco::BeamSpot& s) const {
0151 return vertices(vector<reco::TransientTrack>(), t, s, false, true);
0152 }
0153
0154 vector<TransientVertex> AdaptiveVertexReconstructor::vertices(const vector<reco::TransientTrack>& primaries,
0155 const vector<reco::TransientTrack>& tracks,
0156 const reco::BeamSpot& s) const {
0157 return vertices(primaries, tracks, s, true, true);
0158 }
0159
0160 vector<TransientVertex> AdaptiveVertexReconstructor::vertices(const vector<reco::TransientTrack>& tracks) const {
0161 return vertices(vector<reco::TransientTrack>(), tracks, reco::BeamSpot(), false, false);
0162 }
0163
0164 vector<TransientVertex> AdaptiveVertexReconstructor::vertices(const vector<reco::TransientTrack>& primaries,
0165 const vector<reco::TransientTrack>& tracks,
0166 const reco::BeamSpot& s,
0167 bool has_primaries,
0168 bool usespot) const {
0169 vector<TransientVertex> ret;
0170 set<reco::TransientTrack> remainingtrks;
0171
0172 copy(tracks.begin(), tracks.end(), inserter(remainingtrks, remainingtrks.begin()));
0173
0174 unsigned int n_tracks = remainingtrks.size();
0175
0176
0177 try {
0178 while (remainingtrks.size() > 1) {
0179
0180
0181
0182
0183 const AdaptiveVertexFitter* fitter = theSecondaryFitter;
0184 if (ret.empty()) {
0185 fitter = thePrimaryFitter;
0186 };
0187 vector<reco::TransientTrack> fittrks;
0188 fittrks.reserve(remainingtrks.size());
0189
0190 copy(remainingtrks.begin(), remainingtrks.end(), back_inserter(fittrks));
0191
0192 TransientVertex tmpvtx;
0193 if ((ret.empty()) && has_primaries) {
0194
0195 copy(primaries.begin(), primaries.end(), back_inserter(fittrks));
0196 }
0197 if ((ret.empty()) && usespot) {
0198 tmpvtx = fitter->vertex(fittrks, s);
0199 } else {
0200 tmpvtx = fitter->vertex(fittrks);
0201 }
0202 TransientVertex newvtx = cleanUp(tmpvtx);
0203 ret.push_back(newvtx);
0204 erase(newvtx, remainingtrks, theMinWeight);
0205 if (n_tracks == remainingtrks.size()) {
0206 if (usespot) {
0207
0208 usespot = false;
0209 LogDebug("AdaptiveVertexReconstructor") << "no tracks in vertex. trying again without beamspot constraint!";
0210 continue;
0211 }
0212 LogDebug("AdaptiveVertexReconstructor")
0213 << "all tracks (" << n_tracks << ") would be recycled for next fit. Trying with low threshold!";
0214 erase(newvtx, remainingtrks, 1.e-5);
0215 if (n_tracks == remainingtrks.size()) {
0216 LogDebug("AdaptiveVertexReconstructor") << "low threshold didnt help! "
0217 << "Discontinue procedure!";
0218 break;
0219 }
0220 };
0221
0222
0223 n_tracks = remainingtrks.size();
0224 };
0225 } catch (VertexException& v) {
0226
0227
0228 };
0229
0230 return cleanUpVertices(ret);
0231 }
0232
0233 vector<TransientVertex> AdaptiveVertexReconstructor::cleanUpVertices(const vector<TransientVertex>& old) const {
0234 vector<TransientVertex> ret;
0235 for (vector<TransientVertex>::const_iterator i = old.begin(); i != old.end(); ++i) {
0236 if (!(i->hasTrackWeight())) {
0237 ret.push_back(*i);
0238 continue;
0239 }
0240
0241
0242
0243 int nsig = 0;
0244 TransientVertex::TransientTrackToFloatMap wm = i->weightMap();
0245 for (TransientVertex::TransientTrackToFloatMap::const_iterator w = wm.begin(); w != wm.end(); ++w) {
0246 if (w->second > theWeightThreshold)
0247 nsig++;
0248 }
0249 if (nsig > 1)
0250 ret.push_back(*i);
0251 }
0252
0253 return ret;
0254 }