VertexTools/src/SequentialVertexFitter.cc

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   // FIXME
0011   // hard-coded tracker bounds
0012   // workaround while waiting for Geometry service
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 }  // namespace
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");     //0.01
0068   theMaxStep = thePSet.getParameter<int>("maxNbrOfIterations");  //10
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);  //10
0075   readParameters();
0076 }
0077
0078 template <unsigned int N>
0079 CachingVertex<N> SequentialVertexFitter<N>::vertex(const std::vector<reco::TransientTrack>& tracks) const {
0080   // Linearization Point
0081   GlobalPoint linP = theLinP->getLinearizationPoint(tracks);
0082   if (!insideTrackerBounds(linP))
0083     linP = GlobalPoint(0, 0, 0);
0084
0085   // Initial vertex state, with a very large error matrix
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   // Initial vertex state, with a very small weight matrix
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 // Fit vertex out of a set of RecTracks.
0113 // Uses the specified linearization point.
0114 //
0115 template <unsigned int N>
0116 CachingVertex<N> SequentialVertexFitter<N>::vertex(const std::vector<reco::TransientTrack>& tracks,
0117                                                    const GlobalPoint& linPoint) const {
0118   // Initial vertex state, with a very large error matrix
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 /** Fit vertex out of a set of TransientTracks.
0128    *  The specified BeamSpot will be used as priot, but NOT for the linearization.
0129    * The specified LinearizationPointFinder will be used to find the linearization point.
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     // Linearization Point search if there are more than 1 track
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     // otherwise take the beamspot position.
0149     vtContainer = linearizeTracks(tracks, beamSpotState);
0150   }
0151
0152   return fit(vtContainer, beamSpotState, true);
0153 }
0154
0155 // Fit vertex out of a set of RecTracks.
0156 // Uses the position as both the linearization point AND as prior
0157 // estimate of the vertex position. The error is used for the
0158 // weight of the prior estimate.
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 // Fit vertex out of a set of VertexTracks
0170 // Uses the position and error for the prior estimate of the vertex.
0171 // This position is not used to relinearize the tracks.
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 // Construct a container of VertexTrack from a set of RecTracks.
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 // Construct new a container of VertexTrack with a new linearization point
0198 // and vertex state, from an existing set of VertexTrack, from which only the
0199 // recTracks will be used.
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     //    RefCountedLinearizedTrackState lTrData =
0210     //      theLTrackFactory->linearizedTrackState(linP,
0211     //                  (**i).linearizedTrack()->track());
0212     RefCountedVertexTrack vTrData = theVTrackFactory.vertexTrack(lTrData, state, (**i).weight());
0213     finalTracks.push_back(vTrData);
0214   }
0215   return finalTracks;
0216 }
0217
0218 // The method where the vertex fit is actually done!
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   // main loop through all the VTracks
0237   bool validVertex = true;
0238   int step = 0;
0239   GlobalPoint newPosition = priorVertexPosition;
0240   GlobalPoint previousPosition;
0241   do {
0242     CachingVertex<N> fVertex = initialVertex;
0243     // make new linearized and vertex tracks for the next iteration
0244     if (step != 0)
0245       globalVTracks = reLinearizeTracks(tracks, returnVertex.vertexState());
0246
0247     // update sequentially the vertex estimate
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     // check tracker bounds and NaN in position
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       // reset initial vertex position to (0,0,0) and force new iteration
0270       // if number of steps not exceeded
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>();  // return invalid vertex
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>();  // return invalid vertex
0295   }
0296
0297   // smoothing
0298   returnVertex = theSmoother->smooth(returnVertex);
0299
0300   return returnVertex;
0301 }
0302
0303 template class SequentialVertexFitter<5>;
0304 template class SequentialVertexFitter<6>;