File indexing completed on 2021-02-14 14:28:25
0001 #include "RecoVertex/VertexTools/interface/SequentialVertexFitter.h"
0002 #include "DataFormats/GeometryCommonDetAlgo/interface/GlobalError.h"
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004
0005 #include <algorithm>
0006 using namespace std;
0007 using namespace reco;
0008
0009 namespace {
0010
0011
0012
0013 const float TrackerBoundsRadius = 112;
0014 const float TrackerBoundsHalfLength = 273.5;
0015 bool insideTrackerBounds(const GlobalPoint& point) {
0016 return ((point.transverse() < TrackerBoundsRadius) && (abs(point.z()) < TrackerBoundsHalfLength));
0017 }
0018 }
0019
0020 template <unsigned int N>
0021 SequentialVertexFitter<N>::SequentialVertexFitter(const LinearizationPointFinder& linP,
0022 const VertexUpdator<N>& updator,
0023 const VertexSmoother<N>& smoother,
0024 const AbstractLTSFactory<N>& ltsf)
0025 : theLinP(linP.clone()),
0026 theUpdator(updator.clone()),
0027 theSmoother(smoother.clone()),
0028 theLTrackFactory(ltsf.clone()) {
0029 setDefaultParameters();
0030 }
0031
0032 template <unsigned int N>
0033 SequentialVertexFitter<N>::SequentialVertexFitter(const edm::ParameterSet& pSet,
0034 const LinearizationPointFinder& linP,
0035 const VertexUpdator<N>& updator,
0036 const VertexSmoother<N>& smoother,
0037 const AbstractLTSFactory<N>& ltsf)
0038 : thePSet(pSet),
0039 theLinP(linP.clone()),
0040 theUpdator(updator.clone()),
0041 theSmoother(smoother.clone()),
0042 theLTrackFactory(ltsf.clone()) {
0043 readParameters();
0044 }
0045
0046 template <unsigned int N>
0047 SequentialVertexFitter<N>::SequentialVertexFitter(const SequentialVertexFitter& original) {
0048 thePSet = original.parameterSet();
0049 theLinP = original.linearizationPointFinder()->clone();
0050 theUpdator = original.vertexUpdator()->clone();
0051 theSmoother = original.vertexSmoother()->clone();
0052 theMaxShift = original.maxShift();
0053 theMaxStep = original.maxStep();
0054 theLTrackFactory = original.linearizedTrackStateFactory()->clone();
0055 }
0056
0057 template <unsigned int N>
0058 SequentialVertexFitter<N>::~SequentialVertexFitter() {
0059 delete theLinP;
0060 delete theUpdator;
0061 delete theSmoother;
0062 delete theLTrackFactory;
0063 }
0064
0065 template <unsigned int N>
0066 void SequentialVertexFitter<N>::readParameters() {
0067 theMaxShift = thePSet.getParameter<double>("maxDistance");
0068 theMaxStep = thePSet.getParameter<int>("maxNbrOfIterations");
0069 }
0070
0071 template <unsigned int N>
0072 void SequentialVertexFitter<N>::setDefaultParameters() {
0073 thePSet.addParameter<double>("maxDistance", 0.01);
0074 thePSet.addParameter<int>("maxNbrOfIterations", 10);
0075 readParameters();
0076 }
0077
0078 template <unsigned int N>
0079 CachingVertex<N> SequentialVertexFitter<N>::vertex(const std::vector<reco::TransientTrack>& tracks) const {
0080
0081 GlobalPoint linP = theLinP->getLinearizationPoint(tracks);
0082 if (!insideTrackerBounds(linP))
0083 linP = GlobalPoint(0, 0, 0);
0084
0085
0086 ROOT::Math::SMatrixIdentity id;
0087 AlgebraicSymMatrix33 we(id);
0088 GlobalError error(we * 10000);
0089 VertexState state(linP, error);
0090 std::vector<RefCountedVertexTrack> vtContainer = linearizeTracks(tracks, state);
0091 return fit(vtContainer, state, false);
0092 }
0093
0094 template <unsigned int N>
0095 CachingVertex<N> SequentialVertexFitter<N>::vertex(const std::vector<RefCountedVertexTrack>& tracks,
0096 const reco::BeamSpot& spot) const {
0097 VertexState state(spot);
0098 return fit(tracks, state, true);
0099 }
0100
0101 template <unsigned int N>
0102 CachingVertex<N> SequentialVertexFitter<N>::vertex(const std::vector<RefCountedVertexTrack>& tracks) const {
0103
0104 GlobalPoint linP = tracks[0]->linearizedTrack()->linearizationPoint();
0105 ROOT::Math::SMatrixIdentity id;
0106 AlgebraicSymMatrix33 we(id);
0107 GlobalError error(we * 10000);
0108 VertexState state(linP, error);
0109 return fit(tracks, state, false);
0110 }
0111
0112
0113
0114
0115 template <unsigned int N>
0116 CachingVertex<N> SequentialVertexFitter<N>::vertex(const std::vector<reco::TransientTrack>& tracks,
0117 const GlobalPoint& linPoint) const {
0118
0119 ROOT::Math::SMatrixIdentity id;
0120 AlgebraicSymMatrix33 we(id);
0121 GlobalError error(we * 10000);
0122 VertexState state(linPoint, error);
0123 std::vector<RefCountedVertexTrack> vtContainer = linearizeTracks(tracks, state);
0124 return fit(vtContainer, state, false);
0125 }
0126
0127
0128
0129
0130
0131 template <unsigned int N>
0132 CachingVertex<N> SequentialVertexFitter<N>::vertex(const std::vector<reco::TransientTrack>& tracks,
0133 const BeamSpot& beamSpot) const {
0134 VertexState beamSpotState(beamSpot);
0135 std::vector<RefCountedVertexTrack> vtContainer;
0136
0137 if (tracks.size() > 1) {
0138
0139 GlobalPoint linP = theLinP->getLinearizationPoint(tracks);
0140 if (!insideTrackerBounds(linP))
0141 linP = GlobalPoint(0, 0, 0);
0142 ROOT::Math::SMatrixIdentity id;
0143 AlgebraicSymMatrix33 we(id);
0144 GlobalError error(we * 10000);
0145 VertexState lpState(linP, error);
0146 vtContainer = linearizeTracks(tracks, lpState);
0147 } else {
0148
0149 vtContainer = linearizeTracks(tracks, beamSpotState);
0150 }
0151
0152 return fit(vtContainer, beamSpotState, true);
0153 }
0154
0155
0156
0157
0158
0159
0160 template <unsigned int N>
0161 CachingVertex<N> SequentialVertexFitter<N>::vertex(const std::vector<reco::TransientTrack>& tracks,
0162 const GlobalPoint& priorPos,
0163 const GlobalError& priorError) const {
0164 VertexState state(priorPos, priorError);
0165 std::vector<RefCountedVertexTrack> vtContainer = linearizeTracks(tracks, state);
0166 return fit(vtContainer, state, true);
0167 }
0168
0169
0170
0171
0172
0173 template <unsigned int N>
0174 CachingVertex<N> SequentialVertexFitter<N>::vertex(const std::vector<RefCountedVertexTrack>& tracks,
0175 const GlobalPoint& priorPos,
0176 const GlobalError& priorError) const {
0177 VertexState state(priorPos, priorError);
0178 return fit(tracks, state, true);
0179 }
0180
0181
0182
0183 template <unsigned int N>
0184 vector<typename SequentialVertexFitter<N>::RefCountedVertexTrack> SequentialVertexFitter<N>::linearizeTracks(
0185 const std::vector<reco::TransientTrack>& tracks, const VertexState state) const {
0186 GlobalPoint linP = state.position();
0187 std::vector<RefCountedVertexTrack> finalTracks;
0188 finalTracks.reserve(tracks.size());
0189 for (vector<reco::TransientTrack>::const_iterator i = tracks.begin(); i != tracks.end(); i++) {
0190 RefCountedLinearizedTrackState lTrData = theLTrackFactory->linearizedTrackState(linP, *i);
0191 RefCountedVertexTrack vTrData = theVTrackFactory.vertexTrack(lTrData, state);
0192 finalTracks.push_back(vTrData);
0193 }
0194 return finalTracks;
0195 }
0196
0197
0198
0199
0200
0201 template <unsigned int N>
0202 vector<typename SequentialVertexFitter<N>::RefCountedVertexTrack> SequentialVertexFitter<N>::reLinearizeTracks(
0203 const std::vector<RefCountedVertexTrack>& tracks, const VertexState state) const {
0204 GlobalPoint linP = state.position();
0205 std::vector<RefCountedVertexTrack> finalTracks;
0206 finalTracks.reserve(tracks.size());
0207 for (typename std::vector<RefCountedVertexTrack>::const_iterator i = tracks.begin(); i != tracks.end(); i++) {
0208 RefCountedLinearizedTrackState lTrData = (**i).linearizedTrack()->stateWithNewLinearizationPoint(linP);
0209
0210
0211
0212 RefCountedVertexTrack vTrData = theVTrackFactory.vertexTrack(lTrData, state, (**i).weight());
0213 finalTracks.push_back(vTrData);
0214 }
0215 return finalTracks;
0216 }
0217
0218
0219
0220 template <unsigned int N>
0221 CachingVertex<N> SequentialVertexFitter<N>::fit(const std::vector<RefCountedVertexTrack>& tracks,
0222 const VertexState priorVertex,
0223 bool withPrior) const {
0224 std::vector<RefCountedVertexTrack> initialTracks;
0225 GlobalPoint priorVertexPosition = priorVertex.position();
0226 GlobalError priorVertexError = priorVertex.error();
0227
0228 CachingVertex<N> returnVertex(priorVertexPosition, priorVertexError, initialTracks, 0);
0229 if (withPrior) {
0230 returnVertex = CachingVertex<N>(
0231 priorVertexPosition, priorVertexError, priorVertexPosition, priorVertexError, initialTracks, 0);
0232 }
0233 CachingVertex<N> initialVertex = returnVertex;
0234 std::vector<RefCountedVertexTrack> globalVTracks = tracks;
0235
0236
0237 bool validVertex = true;
0238 int step = 0;
0239 GlobalPoint newPosition = priorVertexPosition;
0240 GlobalPoint previousPosition;
0241 do {
0242 CachingVertex<N> fVertex = initialVertex;
0243
0244 if (step != 0)
0245 globalVTracks = reLinearizeTracks(tracks, returnVertex.vertexState());
0246
0247
0248 for (typename std::vector<RefCountedVertexTrack>::const_iterator i = globalVTracks.begin();
0249 i != globalVTracks.end();
0250 i++) {
0251 fVertex = theUpdator->add(fVertex, *i);
0252 if (!fVertex.isValid())
0253 break;
0254 }
0255
0256 validVertex = fVertex.isValid();
0257
0258 if (validVertex && hasNan(fVertex.position())) {
0259 LogDebug("RecoVertex/SequentialVertexFitter") << "Fitted position is NaN.\n";
0260 validVertex = false;
0261 }
0262
0263 if (validVertex && !insideTrackerBounds(fVertex.position())) {
0264 LogDebug("RecoVertex/SequentialVertexFitter") << "Fitted position is out of tracker bounds.\n";
0265 validVertex = false;
0266 }
0267
0268 if (!validVertex) {
0269
0270
0271 ROOT::Math::SMatrixIdentity id;
0272 AlgebraicSymMatrix33 we(id);
0273 GlobalError error(we * 10000);
0274 fVertex = CachingVertex<N>(GlobalPoint(0, 0, 0), error, initialTracks, 0);
0275 }
0276
0277 previousPosition = newPosition;
0278 newPosition = fVertex.position();
0279
0280 returnVertex = fVertex;
0281 globalVTracks.clear();
0282 step++;
0283 } while ((step != theMaxStep) && (((previousPosition - newPosition).transverse() > theMaxShift) || (!validVertex)));
0284
0285 if (!validVertex) {
0286 LogDebug("RecoVertex/SequentialVertexFitter")
0287 << "Fitted position is invalid (out of tracker bounds or has NaN). Returned vertex is invalid\n";
0288 return CachingVertex<N>();
0289 }
0290
0291 if (step >= theMaxStep) {
0292 LogDebug("RecoVertex/SequentialVertexFitter")
0293 << "The maximum number of steps has been exceeded. Returned vertex is invalid\n";
0294 return CachingVertex<N>();
0295 }
0296
0297
0298 returnVertex = theSmoother->smooth(returnVertex);
0299
0300 return returnVertex;
0301 }
0302
0303 template class SequentialVertexFitter<5>;
0304 template class SequentialVertexFitter<6>;