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 template <unsigned int N>
0010 SequentialVertexFitter<N>::SequentialVertexFitter(const LinearizationPointFinder& linP,
0011                                                   const VertexUpdator<N>& updator,
0012                                                   const VertexSmoother<N>& smoother,
0013                                                   const AbstractLTSFactory<N>& ltsf)
0014     : theLinP(linP.clone()),
0015       theUpdator(updator.clone()),
0016       theSmoother(smoother.clone()),
0017       theLTrackFactory(ltsf.clone()) {
0018   setDefaultParameters();
0019 }
0020
0021 template <unsigned int N>
0022 SequentialVertexFitter<N>::SequentialVertexFitter(const edm::ParameterSet& pSet,
0023                                                   const LinearizationPointFinder& linP,
0024                                                   const VertexUpdator<N>& updator,
0025                                                   const VertexSmoother<N>& smoother,
0026                                                   const AbstractLTSFactory<N>& ltsf)
0027     : thePSet(pSet),
0028       theLinP(linP.clone()),
0029       theUpdator(updator.clone()),
0030       theSmoother(smoother.clone()),
0031       theLTrackFactory(ltsf.clone()) {
0032   readParameters();
0033 }
0034
0035 template <unsigned int N>
0036 SequentialVertexFitter<N>::SequentialVertexFitter(const SequentialVertexFitter& original) {
0037   thePSet = original.parameterSet();
0038   theLinP = original.linearizationPointFinder()->clone();
0039   theUpdator = original.vertexUpdator()->clone();
0040   theSmoother = original.vertexSmoother()->clone();
0041   theMaxShift = original.maxShift();
0042   theMaxStep = original.maxStep();
0043   theLTrackFactory = original.linearizedTrackStateFactory()->clone();
0044 }
0045
0046 template <unsigned int N>
0047 SequentialVertexFitter<N>::~SequentialVertexFitter() {
0048   delete theLinP;
0049   delete theUpdator;
0050   delete theSmoother;
0051   delete theLTrackFactory;
0052 }
0053
0054 template <unsigned int N>
0055 void SequentialVertexFitter<N>::readParameters() {
0056   theMaxShift = thePSet.getParameter<double>("maxDistance");     //0.01
0057   theMaxStep = thePSet.getParameter<int>("maxNbrOfIterations");  //10
0058 }
0059
0060 template <unsigned int N>
0061 void SequentialVertexFitter<N>::setDefaultParameters() {
0062   thePSet.addParameter<double>("maxDistance", 0.01);
0063   thePSet.addParameter<int>("maxNbrOfIterations", 10);  //10
0064   readParameters();
0065 }
0066
0067 template <unsigned int N>
0068 CachingVertex<N> SequentialVertexFitter<N>::vertex(const std::vector<reco::TransientTrack>& tracks) const {
0069   // Linearization Point
0070   GlobalPoint linP = theLinP->getLinearizationPoint(tracks);
0071   if (!insideTrackerBounds(linP))
0072     linP = GlobalPoint(0, 0, 0);
0073
0074   // Initial vertex state, with a very large error matrix
0075   ROOT::Math::SMatrixIdentity id;
0076   AlgebraicSymMatrix33 we(id);
0077   GlobalError error(we * 10000);
0078   VertexState state(linP, error);
0079   std::vector<RefCountedVertexTrack> vtContainer = linearizeTracks(tracks, state);
0080   return fit(vtContainer, state, false);
0081 }
0082
0083 template <unsigned int N>
0084 CachingVertex<N> SequentialVertexFitter<N>::vertex(const std::vector<RefCountedVertexTrack>& tracks,
0085                                                    const reco::BeamSpot& spot) const {
0086   VertexState state(spot);
0087   return fit(tracks, state, true);
0088 }
0089
0090 template <unsigned int N>
0091 CachingVertex<N> SequentialVertexFitter<N>::vertex(const std::vector<RefCountedVertexTrack>& tracks) const {
0092   // Initial vertex state, with a very small weight matrix
0093   GlobalPoint linP = tracks[0]->linearizedTrack()->linearizationPoint();
0094   ROOT::Math::SMatrixIdentity id;
0095   AlgebraicSymMatrix33 we(id);
0096   GlobalError error(we * 10000);
0097   VertexState state(linP, error);
0098   return fit(tracks, state, false);
0099 }
0100
0101 // Fit vertex out of a set of RecTracks.
0102 // Uses the specified linearization point.
0103 //
0104 template <unsigned int N>
0105 CachingVertex<N> SequentialVertexFitter<N>::vertex(const std::vector<reco::TransientTrack>& tracks,
0106                                                    const GlobalPoint& linPoint) const {
0107   // Initial vertex state, with a very large error matrix
0108   ROOT::Math::SMatrixIdentity id;
0109   AlgebraicSymMatrix33 we(id);
0110   GlobalError error(we * 10000);
0111   VertexState state(linPoint, error);
0112   std::vector<RefCountedVertexTrack> vtContainer = linearizeTracks(tracks, state);
0113   return fit(vtContainer, state, false);
0114 }
0115
0116 /** Fit vertex out of a set of TransientTracks.
0117    *  The specified BeamSpot will be used as priot, but NOT for the linearization.
0118    * The specified LinearizationPointFinder will be used to find the linearization point.
0119    */
0120 template <unsigned int N>
0121 CachingVertex<N> SequentialVertexFitter<N>::vertex(const std::vector<reco::TransientTrack>& tracks,
0122                                                    const BeamSpot& beamSpot) const {
0123   VertexState beamSpotState(beamSpot);
0124   std::vector<RefCountedVertexTrack> vtContainer;
0125
0126   if (tracks.size() > 1) {
0127     // Linearization Point search if there are more than 1 track
0128     GlobalPoint linP = theLinP->getLinearizationPoint(tracks);
0129     if (!insideTrackerBounds(linP))
0130       linP = GlobalPoint(0, 0, 0);
0131     ROOT::Math::SMatrixIdentity id;
0132     AlgebraicSymMatrix33 we(id);
0133     GlobalError error(we * 10000);
0134     VertexState lpState(linP, error);
0135     vtContainer = linearizeTracks(tracks, lpState);
0136   } else {
0137     // otherwise take the beamspot position.
0138     vtContainer = linearizeTracks(tracks, beamSpotState);
0139   }
0140
0141   return fit(vtContainer, beamSpotState, true);
0142 }
0143
0144 // Fit vertex out of a set of RecTracks.
0145 // Uses the position as both the linearization point AND as prior
0146 // estimate of the vertex position. The error is used for the
0147 // weight of the prior estimate.
0148 //
0149 template <unsigned int N>
0150 CachingVertex<N> SequentialVertexFitter<N>::vertex(const std::vector<reco::TransientTrack>& tracks,
0151                                                    const GlobalPoint& priorPos,
0152                                                    const GlobalError& priorError) const {
0153   VertexState state(priorPos, priorError);
0154   std::vector<RefCountedVertexTrack> vtContainer = linearizeTracks(tracks, state);
0155   return fit(vtContainer, state, true);
0156 }
0157
0158 // Fit vertex out of a set of VertexTracks
0159 // Uses the position and error for the prior estimate of the vertex.
0160 // This position is not used to relinearize the tracks.
0161 //
0162 template <unsigned int N>
0163 CachingVertex<N> SequentialVertexFitter<N>::vertex(const std::vector<RefCountedVertexTrack>& tracks,
0164                                                    const GlobalPoint& priorPos,
0165                                                    const GlobalError& priorError) const {
0166   VertexState state(priorPos, priorError);
0167   return fit(tracks, state, true);
0168 }
0169
0170 // Construct a container of VertexTrack from a set of RecTracks.
0171 //
0172 template <unsigned int N>
0173 vector<typename SequentialVertexFitter<N>::RefCountedVertexTrack> SequentialVertexFitter<N>::linearizeTracks(
0174     const std::vector<reco::TransientTrack>& tracks, const VertexState state) const {
0175   GlobalPoint linP = state.position();
0176   std::vector<RefCountedVertexTrack> finalTracks;
0177   finalTracks.reserve(tracks.size());
0178   for (vector<reco::TransientTrack>::const_iterator i = tracks.begin(); i != tracks.end(); i++) {
0179     RefCountedLinearizedTrackState lTrData = theLTrackFactory->linearizedTrackState(linP, *i);
0180     RefCountedVertexTrack vTrData = theVTrackFactory.vertexTrack(lTrData, state);
0181     finalTracks.push_back(vTrData);
0182   }
0183   return finalTracks;
0184 }
0185
0186 // Construct new a container of VertexTrack with a new linearization point
0187 // and vertex state, from an existing set of VertexTrack, from which only the
0188 // recTracks will be used.
0189 //
0190 template <unsigned int N>
0191 vector<typename SequentialVertexFitter<N>::RefCountedVertexTrack> SequentialVertexFitter<N>::reLinearizeTracks(
0192     const std::vector<RefCountedVertexTrack>& tracks, const VertexState state) const {
0193   GlobalPoint linP = state.position();
0194   std::vector<RefCountedVertexTrack> finalTracks;
0195   finalTracks.reserve(tracks.size());
0196   for (typename std::vector<RefCountedVertexTrack>::const_iterator i = tracks.begin(); i != tracks.end(); i++) {
0197     RefCountedLinearizedTrackState lTrData = (**i).linearizedTrack()->stateWithNewLinearizationPoint(linP);
0198     //    RefCountedLinearizedTrackState lTrData =
0199     //      theLTrackFactory->linearizedTrackState(linP,
0200     //                  (**i).linearizedTrack()->track());
0201     RefCountedVertexTrack vTrData = theVTrackFactory.vertexTrack(lTrData, state, (**i).weight());
0202     finalTracks.push_back(vTrData);
0203   }
0204   return finalTracks;
0205 }
0206
0207 // The method where the vertex fit is actually done!
0208 //
0209 template <unsigned int N>
0210 CachingVertex<N> SequentialVertexFitter<N>::fit(const std::vector<RefCountedVertexTrack>& tracks,
0211                                                 const VertexState priorVertex,
0212                                                 bool withPrior) const {
0213   std::vector<RefCountedVertexTrack> initialTracks;
0214   GlobalPoint priorVertexPosition = priorVertex.position();
0215   GlobalError priorVertexError = priorVertex.error();
0216
0217   CachingVertex<N> returnVertex(priorVertexPosition, priorVertexError, initialTracks, 0);
0218   if (withPrior) {
0219     returnVertex = CachingVertex<N>(
0220         priorVertexPosition, priorVertexError, priorVertexPosition, priorVertexError, initialTracks, 0);
0221   }
0222   CachingVertex<N> initialVertex = returnVertex;
0223   std::vector<RefCountedVertexTrack> globalVTracks = tracks;
0224
0225   // main loop through all the VTracks
0226   bool validVertex = true;
0227   int step = 0;
0228   GlobalPoint newPosition = priorVertexPosition;
0229   GlobalPoint previousPosition;
0230   do {
0231     CachingVertex<N> fVertex = initialVertex;
0232     // make new linearized and vertex tracks for the next iteration
0233     if (step != 0)
0234       globalVTracks = reLinearizeTracks(tracks, returnVertex.vertexState());
0235
0236     // update sequentially the vertex estimate
0237     for (typename std::vector<RefCountedVertexTrack>::const_iterator i = globalVTracks.begin();
0238          i != globalVTracks.end();
0239          i++) {
0240       fVertex = theUpdator->add(fVertex, *i);
0241       if (!fVertex.isValid())
0242         break;
0243     }
0244
0245     validVertex = fVertex.isValid();
0246     // check tracker bounds and NaN in position
0247     if (validVertex && hasNan(fVertex.position())) {
0248       LogDebug("RecoVertex/SequentialVertexFitter") << "Fitted position is NaN.\n";
0249       validVertex = false;
0250     }
0251
0252     if (validVertex && !insideTrackerBounds(fVertex.position())) {
0253       LogDebug("RecoVertex/SequentialVertexFitter") << "Fitted position is out of tracker bounds.\n";
0254       validVertex = false;
0255     }
0256
0257     if (!validVertex) {
0258       // reset initial vertex position to (0,0,0) and force new iteration
0259       // if number of steps not exceeded
0260       ROOT::Math::SMatrixIdentity id;
0261       AlgebraicSymMatrix33 we(id);
0262       GlobalError error(we * 10000);
0263       fVertex = CachingVertex<N>(GlobalPoint(0, 0, 0), error, initialTracks, 0);
0264     }
0265
0266     previousPosition = newPosition;
0267     newPosition = fVertex.position();
0268
0269     returnVertex = fVertex;
0270     globalVTracks.clear();
0271     step++;
0272   } while ((step != theMaxStep) && (((previousPosition - newPosition).transverse() > theMaxShift) || (!validVertex)));
0273
0274   if (!validVertex) {
0275     LogDebug("RecoVertex/SequentialVertexFitter")
0276         << "Fitted position is invalid (out of tracker bounds or has NaN). Returned vertex is invalid\n";
0277     return CachingVertex<N>();  // return invalid vertex
0278   }
0279
0280   if (step >= theMaxStep) {
0281     LogDebug("RecoVertex/SequentialVertexFitter")
0282         << "The maximum number of steps has been exceeded. Returned vertex is invalid\n";
0283     return CachingVertex<N>();  // return invalid vertex
0284   }
0285
0286   // smoothing
0287   returnVertex = theSmoother->smooth(returnVertex);
0288
0289   return returnVertex;
0290 }
0291
0292 template class SequentialVertexFitter<5>;
0293 template class SequentialVertexFitter<6>;