GaussianSumVertexFit/src/GsfVertexSmoother.cc

0001 #include "RecoVertex/GaussianSumVertexFit/interface/GsfVertexSmoother.h"
0002 #include "RecoVertex/GaussianSumVertexFit/interface/BasicMultiVertexState.h"
0003 #include "RecoVertex/GaussianSumVertexFit/interface/MultiRefittedTS.h"
0004 #include "RecoVertex/VertexPrimitives/interface/VertexException.h"
0005
0006 GsfVertexSmoother::GsfVertexSmoother(bool limit, const GsfVertexMerger* merger) : limitComponents(limit) {
0007   if (limitComponents)
0008     theMerger = merger->clone();
0009 }
0010
0011 CachingVertex<5> GsfVertexSmoother::smooth(const CachingVertex<5>& vertex) const {
0012   std::vector<RefCountedVertexTrack> tracks = vertex.tracks();
0013   int numberOfTracks = tracks.size();
0014   if (numberOfTracks < 1)
0015     return vertex;
0016
0017   // Initial vertex for ascending fit
0018   GlobalPoint priorVertexPosition = tracks[0]->linearizedTrack()->linearizationPoint();
0019   AlgebraicSymMatrix33 we = ROOT::Math::SMatrixIdentity();
0020   GlobalError priorVertexError(we * 10000);
0021
0022   std::vector<RefCountedVertexTrack> initialTracks;
0023   CachingVertex<5> fitVertex(priorVertexPosition, priorVertexError, initialTracks, 0);
0024   //In case prior vertex was used.
0025   if (vertex.hasPrior()) {
0026     const VertexState& priorVertexState = vertex.priorVertexState();
0027     fitVertex = CachingVertex<5>(priorVertexState, priorVertexState, initialTracks, 0);
0028   }
0029
0030   // vertices from ascending fit
0031   std::vector<CachingVertex<5> > ascendingFittedVertices;
0032   ascendingFittedVertices.reserve(numberOfTracks);
0033   ascendingFittedVertices.push_back(fitVertex);
0034
0035   // ascending fit
0036   for (std::vector<RefCountedVertexTrack>::const_iterator i = tracks.begin(); i != (tracks.end() - 1); ++i) {
0037     fitVertex = theUpdator.add(fitVertex, *i);
0038     if (limitComponents)
0039       fitVertex = theMerger->merge(fitVertex);
0040     ascendingFittedVertices.push_back(fitVertex);
0041   }
0042
0043   // Initial vertex for descending fit
0044   priorVertexPosition = tracks[0]->linearizedTrack()->linearizationPoint();
0045   priorVertexError = GlobalError(we * 10000);
0046   fitVertex = CachingVertex<5>(priorVertexPosition, priorVertexError, initialTracks, 0);
0047
0048   // vertices from descending fit
0049   std::vector<CachingVertex<5> > descendingFittedVertices;
0050   descendingFittedVertices.reserve(numberOfTracks);
0051   descendingFittedVertices.push_back(fitVertex);
0052
0053   // descending fit
0054   for (std::vector<RefCountedVertexTrack>::const_iterator i = (tracks.end() - 1); i != (tracks.begin()); --i) {
0055     fitVertex = theUpdator.add(fitVertex, *i);
0056     if (limitComponents)
0057       fitVertex = theMerger->merge(fitVertex);
0058     descendingFittedVertices.insert(descendingFittedVertices.begin(), fitVertex);
0059   }
0060
0061   std::vector<RefCountedVertexTrack> newTracks;
0062   double smoothedChi2 = 0.;  // Smoothed chi2
0063
0064   // Track refit
0065   for (std::vector<RefCountedVertexTrack>::const_iterator i = tracks.begin(); i != tracks.end(); i++) {
0066     int indexNumber = i - tracks.begin();
0067     //First, combine the vertices:
0068     VertexState meanedVertex = meanVertex(ascendingFittedVertices[indexNumber].vertexState(),
0069                                           descendingFittedVertices[indexNumber].vertexState());
0070     if (limitComponents)
0071       meanedVertex = theMerger->merge(meanedVertex);
0072     // Add the vertex and smooth the track:
0073     TrackChi2Pair thePair = vertexAndTrackUpdate(meanedVertex, *i, vertex.position());
0074     smoothedChi2 += thePair.second.second;
0075     newTracks.push_back(theVTFactory.vertexTrack((**i).linearizedTrack(),
0076                                                  vertex.vertexState(),
0077                                                  thePair.first,
0078                                                  thePair.second.second,
0079                                                  AlgebraicSymMatrixOO(),
0080                                                  (**i).weight()));
0081   }
0082
0083   if (vertex.hasPrior()) {
0084     smoothedChi2 += priorVertexChi2(vertex.priorVertexState(), vertex.vertexState());
0085     return CachingVertex<5>(vertex.priorVertexState(), vertex.vertexState(), newTracks, smoothedChi2);
0086   } else {
0087     return CachingVertex<5>(vertex.vertexState(), newTracks, smoothedChi2);
0088   }
0089 }
0090
0091 GsfVertexSmoother::TrackChi2Pair GsfVertexSmoother::vertexAndTrackUpdate(const VertexState& oldVertex,
0092                                                                          const RefCountedVertexTrack track,
0093                                                                          const GlobalPoint& referencePosition) const {
0094   VSC prevVtxComponents = oldVertex.components();
0095
0096   if (prevVtxComponents.empty()) {
0097     throw VertexException("GsfVertexSmoother::(Previous) Vertex to update has no components");
0098   }
0099
0100   LTC ltComponents = track->linearizedTrack()->components();
0101   if (ltComponents.empty()) {
0102     throw VertexException("GsfVertexSmoother::Track to add to vertex has no components");
0103   }
0104   float trackWeight = track->weight();
0105
0106   std::vector<RefittedTrackComponent> newTrackComponents;
0107   newTrackComponents.reserve(prevVtxComponents.size() * ltComponents.size());
0108
0109   for (VSC::iterator vertexCompIter = prevVtxComponents.begin(); vertexCompIter != prevVtxComponents.end();
0110        vertexCompIter++) {
0111     for (LTC::iterator trackCompIter = ltComponents.begin(); trackCompIter != ltComponents.end(); trackCompIter++) {
0112       newTrackComponents.push_back(createNewComponent(*vertexCompIter, *trackCompIter, trackWeight));
0113     }
0114   }
0115
0116   return assembleTrackComponents(newTrackComponents, referencePosition);
0117 }
0118
0119 /**
0120  * This method assembles all the components of the refitted track into one refitted track state,
0121  * normalizing the components. Also, it adds the chi2 track-components increments.
0122  */
0123
0124 GsfVertexSmoother::TrackChi2Pair GsfVertexSmoother::assembleTrackComponents(
0125     const std::vector<GsfVertexSmoother::RefittedTrackComponent>& trackComponents,
0126     const GlobalPoint& referencePosition) const {
0127   //renormalize weights
0128
0129   double totalWeight = 0.;
0130   double totalVtxChi2 = 0., totalTrkChi2 = 0.;
0131
0132   for (std::vector<RefittedTrackComponent>::const_iterator iter = trackComponents.begin();
0133        iter != trackComponents.end();
0134        ++iter) {
0135     totalWeight += iter->first.second;
0136     totalVtxChi2 += iter->second.first * iter->first.second;
0137     totalTrkChi2 += iter->second.second * iter->first.second;
0138   }
0139
0140   totalVtxChi2 /= totalWeight;
0141   totalTrkChi2 /= totalWeight;
0142
0143   std::vector<RefCountedRefittedTrackState> reWeightedRTSC;
0144   reWeightedRTSC.reserve(trackComponents.size());
0145
0146   for (std::vector<RefittedTrackComponent>::const_iterator iter = trackComponents.begin();
0147        iter != trackComponents.end();
0148        ++iter) {
0149     if (iter->second.first != 0) {
0150       reWeightedRTSC.push_back(iter->first.first->stateWithNewWeight(iter->second.first / totalWeight));
0151     }
0152   }
0153
0154   RefCountedRefittedTrackState finalRTS =
0155       RefCountedRefittedTrackState(new MultiRefittedTS(reWeightedRTSC, referencePosition));
0156   return TrackChi2Pair(finalRTS, VtxTrkChi2Pair(totalVtxChi2, totalTrkChi2));
0157 }
0158
0159 /**
0160    * This method does the smoothing of one track component with one vertex component.
0161    * And the track-component-chi2 increment and weight of new component in mixture.
0162    */
0163
0164 GsfVertexSmoother::RefittedTrackComponent GsfVertexSmoother::createNewComponent(
0165     const VertexState& oldVertex, const RefCountedLinearizedTrackState linTrack, float trackWeight) const {
0166   int sign = +1;
0167
0168   // Weight of the component in the mixture (non-normalized)
0169   double weightInMixture = theWeightCalculator.calculate(oldVertex, linTrack, 1000000000.);
0170
0171   // position estimate of the component
0172   VertexState newVertex = kalmanVertexUpdator.positionUpdate(oldVertex, linTrack, trackWeight, sign);
0173
0174   KalmanVertexTrackUpdator<5>::trackMatrixPair thePair =
0175       theVertexTrackUpdator.trackRefit(newVertex, linTrack, trackWeight);
0176
0177   //Chi**2 contribution of the track component
0178   double vtxChi2 = helper.vertexChi2(oldVertex, newVertex);
0179   std::pair<bool, double> trkCi2 = helper.trackParameterChi2(linTrack, thePair.first);
0180
0181   return RefittedTrackComponent(TrackWeightPair(thePair.first, weightInMixture),
0182                                 VtxTrkChi2Pair(vtxChi2, trkCi2.second));
0183 }
0184
0185 VertexState GsfVertexSmoother::meanVertex(const VertexState& vertexA, const VertexState& vertexB) const {
0186   std::vector<VertexState> vsCompA = vertexA.components();
0187   std::vector<VertexState> vsCompB = vertexB.components();
0188   std::vector<VertexState> finalVS;
0189   finalVS.reserve(vsCompA.size() * vsCompB.size());
0190   for (std::vector<VertexState>::iterator iA = vsCompA.begin(); iA != vsCompA.end(); ++iA) {
0191     for (std::vector<VertexState>::iterator iB = vsCompB.begin(); iB != vsCompB.end(); ++iB) {
0192       AlgebraicSymMatrix33 newWeight = iA->weight().matrix() + iB->weight().matrix();
0193       AlgebraicVector3 newWtP = iA->weightTimesPosition() + iB->weightTimesPosition();
0194       double newWeightInMixture = iA->weightInMixture() * iB->weightInMixture();
0195       finalVS.push_back(VertexState(newWtP, newWeight, newWeightInMixture));
0196     }
0197   }
0198 #ifndef CMS_NO_COMPLEX_RETURNS
0199   return VertexState(new BasicMultiVertexState(finalVS));
0200 #else
0201   VertexState theFinalVM(new BasicMultiVertexState(finalVS));
0202   return theFinalVM;
0203 #endif
0204 }
0205
0206 double GsfVertexSmoother::priorVertexChi2(const VertexState priorVertex, const VertexState fittedVertex) const {
0207   std::vector<VertexState> priorVertexComp = priorVertex.components();
0208   std::vector<VertexState> fittedVertexComp = fittedVertex.components();
0209   double vetexChi2 = 0.;
0210   for (std::vector<VertexState>::iterator pvI = priorVertexComp.begin(); pvI != priorVertexComp.end(); ++pvI) {
0211     for (std::vector<VertexState>::iterator fvI = fittedVertexComp.begin(); fvI != fittedVertexComp.end(); ++fvI) {
0212       vetexChi2 += (pvI->weightInMixture()) * (fvI->weightInMixture()) * helper.vertexChi2(*pvI, *fvI);
0213     }
0214   }
0215   return vetexChi2;
0216 }