GaussianSumVertexFit/src/GsfVertexUpdator.cc

0001 #include "RecoVertex/GaussianSumVertexFit/interface/GsfVertexUpdator.h"
0002 #include "RecoVertex/GaussianSumVertexFit/interface/BasicMultiVertexState.h"
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004 #include <cfloat>
0005
0006 GsfVertexUpdator::GsfVertexUpdator(bool limit, const GsfVertexMerger* merger) : limitComponents(limit) {
0007   if (limitComponents)
0008     theMerger = merger->clone();
0009 }
0010
0011 CachingVertex<5> GsfVertexUpdator::add(const CachingVertex<5>& oldVertex, const RefCountedVertexTrack track) const {
0012   VSC prevVtxComponents = oldVertex.vertexState().components();
0013
0014   //   cout << "GsfVertexUpdator::Add new Track with "
0015   //       << track->linearizedTrack()->components().size() << " components to vertex of "
0016   //       << prevVtxComponents.size() << " components.\n";
0017
0018   if (prevVtxComponents.empty()) {
0019     throw VertexException("GsfVertexUpdator::(Previous) Vertex to update has no components");
0020   }
0021
0022   LTC ltComponents = track->linearizedTrack()->components();
0023   if (ltComponents.empty()) {
0024     throw VertexException("GsfVertexUpdator::Track to add to vertex has no components");
0025   }
0026
0027   if ((ltComponents.size() == 1) && (prevVtxComponents.size() == 1))
0028     return kalmanVertexUpdator.add(oldVertex, track);
0029
0030   float trackWeight = track->weight();
0031
0032   std::vector<VertexComponent> newVertexComponents;
0033   newVertexComponents.reserve(prevVtxComponents.size() * ltComponents.size());
0034
0035   //     for (LTC::iterator trackCompIter = ltComponents.begin();
0036   //    trackCompIter != ltComponents.end(); trackCompIter++ ) {
0037   //   cout <<(**trackCompIter).state().globalPosition()<<endl;
0038   //     }
0039
0040   for (VSC::iterator vertexCompIter = prevVtxComponents.begin(); vertexCompIter != prevVtxComponents.end();
0041        vertexCompIter++) {
0042     for (LTC::iterator trackCompIter = ltComponents.begin(); trackCompIter != ltComponents.end(); trackCompIter++) {
0043       newVertexComponents.push_back(createNewComponent(*vertexCompIter, *trackCompIter, trackWeight, +1));
0044       // return invalid vertex in case one of the updated vertex-components is invalid
0045       if (!newVertexComponents.back().first.isValid())
0046         return CachingVertex<5>();
0047     }
0048   }
0049   //   cout << "updator components: "<<newVertexComponents.size()<<endl;
0050
0051   // Update tracks vector
0052
0053   std::vector<RefCountedVertexTrack> newVertexTracks = oldVertex.tracks();
0054   newVertexTracks.push_back(track);
0055   //   cout << "a \n ";
0056
0057   // Assemble VertexStates and compute Chi**2
0058
0059   VertexChi2Pair vertexChi2Pair = assembleVertexComponents(newVertexComponents);
0060   //   cout << "b \n ";
0061   VertexState newVertexState = vertexChi2Pair.first;
0062   //   cout << "c \n ";
0063   double chi2 = oldVertex.totalChiSquared() + vertexChi2Pair.second;
0064
0065   // Merge:
0066   if (limitComponents)
0067     newVertexState = theMerger->merge(newVertexState);
0068
0069   if (oldVertex.hasPrior()) {
0070     return CachingVertex<5>(oldVertex.priorPosition(),
0071                             oldVertex.priorError(),
0072                             newVertexState.weightTimesPosition(),
0073                             newVertexState.weight(),
0074                             newVertexTracks,
0075                             chi2);
0076   } else {
0077     return CachingVertex<5>(newVertexState, newVertexTracks, chi2);
0078   }
0079 }
0080
0081 CachingVertex<5> GsfVertexUpdator::remove(const CachingVertex<5>& oldVertex, const RefCountedVertexTrack track) const {
0082   throw VertexException("GsfVertexUpdator::Remove Methode not yet done");
0083   //  return CachingVertex<5>();
0084 }
0085
0086 /**
0087    * Where one component of the previous vertex gets updated with one component of
0088    *  the track.
0089    */
0090
0091 GsfVertexUpdator::VertexComponent GsfVertexUpdator::createNewComponent(const VertexState& oldVertex,
0092                                                                        const RefCountedLinearizedTrackState linTrack,
0093                                                                        float weight,
0094                                                                        int sign) const {
0095   if (abs(sign) != 1)
0096     throw VertexException("GsfVertexUpdator::sign not equal to 1.");
0097
0098   if (sign == -1)
0099     throw VertexException("GsfVertexUpdator::sign of -1 not yet implemented.");
0100
0101   // Weight of the component in the mixture (non-normalized)
0102   double weightInMixture = theWeightCalculator.calculate(oldVertex, linTrack, 1.E9);
0103   if (weightInMixture < 0.)
0104     return VertexComponent(VertexState(), WeightChi2Pair(0., 0.));
0105
0106   // position estimate of the component
0107   VertexState newVertex = kalmanVertexUpdator.positionUpdate(oldVertex, linTrack, weight, sign);
0108   if (!newVertex.isValid())
0109     return VertexComponent(newVertex, WeightChi2Pair(0., 0.));
0110
0111   //Chi**2 contribution of the component
0112   std::pair<bool, double> chi2P = kalmanVertexUpdator.chi2Increment(oldVertex, newVertex, linTrack, weight);
0113   if (!chi2P.first)
0114     return VertexComponent(VertexState(), WeightChi2Pair(0., 0.));
0115   //         cout << "Update: "<<oldVertex.position()<<" "<<newVertex.position()<<" "<<chi2
0116   //         <<" "<<linTrack->weightInMixture()<<" "<<weightInMixture<<endl;
0117
0118   return VertexComponent(newVertex, WeightChi2Pair(weightInMixture, chi2P.second));
0119 }
0120
0121 GsfVertexUpdator::VertexChi2Pair GsfVertexUpdator::assembleVertexComponents(
0122     const std::vector<GsfVertexUpdator::VertexComponent>& newVertexComponents) const {
0123   VSC vertexComponents;
0124   vertexComponents.reserve(newVertexComponents.size());
0125
0126   //renormalize weights
0127   // cout << "assemble "<<newVertexComponents.size()<<endl;
0128   double totalWeight = 0.;
0129   double totalChi2 = 0.;
0130
0131   for (std::vector<VertexComponent>::const_iterator iter = newVertexComponents.begin();
0132        iter != newVertexComponents.end();
0133        iter++) {
0134     totalWeight += iter->second.first;
0135   }
0136   // cout << "totalWeight "<<totalWeight<<endl;
0137   if (totalWeight < DBL_MIN) {
0138     edm::LogWarning("GsfVertexUpdator") << "Updated Vertex has total weight of 0. "
0139                                         << "The track is probably very far away.";
0140     return VertexChi2Pair(VertexState(), 0.);
0141   }
0142
0143   for (std::vector<VertexComponent>::const_iterator iter = newVertexComponents.begin();
0144        iter != newVertexComponents.end();
0145        iter++) {
0146     double weight = iter->second.first / totalWeight;
0147     if (iter->second.first > DBL_MIN) {
0148       vertexComponents.push_back(VertexState(iter->first.weightTimesPosition(), iter->first.weight(), weight));
0149       totalChi2 += iter->second.second * weight;
0150     }
0151   }
0152   // cout << "totalChi2 "<<totalChi2<<endl;
0153   // cout << "vertexComponents "<<vertexComponents.size()<<endl;
0154
0155   if (vertexComponents.empty()) {
0156     edm::LogWarning("GsfVertexUpdator") << "No Vertex State left after reweighting.";
0157     return VertexChi2Pair(VertexState(), 0.);
0158   }
0159
0160   return VertexChi2Pair(VertexState(new BasicMultiVertexState(vertexComponents)), totalChi2);
0161 }