File indexing completed on 2023-10-25 10:03:18
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
0015
0016
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
0036
0037
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
0045 if (!newVertexComponents.back().first.isValid())
0046 return CachingVertex<5>();
0047 }
0048 }
0049
0050
0051
0052
0053 std::vector<RefCountedVertexTrack> newVertexTracks = oldVertex.tracks();
0054 newVertexTracks.push_back(track);
0055
0056
0057
0058
0059 VertexChi2Pair vertexChi2Pair = assembleVertexComponents(newVertexComponents);
0060
0061 VertexState newVertexState = vertexChi2Pair.first;
0062
0063 double chi2 = oldVertex.totalChiSquared() + vertexChi2Pair.second;
0064
0065
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
0084 }
0085
0086
0087
0088
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
0102 double weightInMixture = theWeightCalculator.calculate(oldVertex, linTrack, 1.E9);
0103 if (weightInMixture < 0.)
0104 return VertexComponent(VertexState(), WeightChi2Pair(0., 0.));
0105
0106
0107 VertexState newVertex = kalmanVertexUpdator.positionUpdate(oldVertex, linTrack, weight, sign);
0108 if (!newVertex.isValid())
0109 return VertexComponent(newVertex, WeightChi2Pair(0., 0.));
0110
0111
0112 std::pair<bool, double> chi2P = kalmanVertexUpdator.chi2Increment(oldVertex, newVertex, linTrack, weight);
0113 if (!chi2P.first)
0114 return VertexComponent(VertexState(), WeightChi2Pair(0., 0.));
0115
0116
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
0127
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
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
0153
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 }