File indexing completed on 2024-04-06 12:31:32
0001 #include "TrackingTools/GsfTracking/interface/PosteriorWeightsCalculator.h"
0002
0003 #include "TrackingTools/PatternTools/interface/MeasurementExtractor.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005 #include "DataFormats/TrackingRecHit/interface/KfComponentsHolder.h"
0006 #include "DataFormats/Math/interface/invertPosDefMatrix.h"
0007 #include "DataFormats/Math/interface/ProjectMatrix.h"
0008
0009 #include <cfloat>
0010
0011 std::vector<double> PosteriorWeightsCalculator::weights(const TrackingRecHit& recHit) const {
0012 switch (recHit.dimension()) {
0013 case 1:
0014 return weights<1>(recHit);
0015 case 2:
0016 return weights<2>(recHit);
0017 case 3:
0018 return weights<3>(recHit);
0019 case 4:
0020 return weights<4>(recHit);
0021 case 5:
0022 return weights<5>(recHit);
0023 }
0024 throw cms::Exception("Error: rechit of size not 1,2,3,4,5");
0025 }
0026
0027 template <unsigned int D>
0028 std::vector<double> PosteriorWeightsCalculator::weights(const TrackingRecHit& recHit) const {
0029 typedef typename AlgebraicROOTObject<D, 5>::Matrix MatD5;
0030 typedef typename AlgebraicROOTObject<5, D>::Matrix Mat5D;
0031 typedef typename AlgebraicROOTObject<D, D>::SymMatrix SMatDD;
0032 typedef typename AlgebraicROOTObject<D>::Vector VecD;
0033 using ROOT::Math::SMatrixNoInit;
0034
0035 std::vector<double> weights;
0036 if (predictedComponents.empty()) {
0037 edm::LogError("EmptyPredictedComponents") << "a multi state is empty. cannot compute any weight.";
0038 return weights;
0039 }
0040 weights.reserve(predictedComponents.size());
0041
0042 std::vector<double> detRs;
0043 detRs.reserve(predictedComponents.size());
0044 std::vector<double> chi2s;
0045 chi2s.reserve(predictedComponents.size());
0046
0047 VecD r, rMeas;
0048 SMatDD V(SMatrixNoInit{}), R(SMatrixNoInit{});
0049 ProjectMatrix<double, 5, D> p;
0050
0051
0052
0053
0054 double chi2Min(DBL_MAX);
0055 for (unsigned int i = 0; i < predictedComponents.size(); i++) {
0056 KfComponentsHolder holder;
0057 auto const& x = predictedComponents[i].localParameters().vector();
0058 holder.template setup<D>(&r, &V, &p, &rMeas, &R, x, predictedComponents[i].localError().matrix());
0059 recHit.getKfComponents(holder);
0060
0061 r -= rMeas;
0062 R += V;
0063
0064 double detR;
0065 if (!R.Det2(detR)) {
0066 edm::LogError("PosteriorWeightsCalculator") << "PosteriorWeightsCalculator: determinant failed";
0067 return std::vector<double>();
0068 }
0069 detRs.push_back(detR);
0070
0071 bool ok = invertPosDefMatrix(R);
0072 if (!ok) {
0073 edm::LogError("PosteriorWeightsCalculator") << "PosteriorWeightsCalculator: inversion failed";
0074 return std::vector<double>();
0075 }
0076 double chi2 = ROOT::Math::Similarity(r, R);
0077 chi2s.push_back(chi2);
0078 if (chi2 < chi2Min)
0079 chi2Min = chi2;
0080 }
0081
0082 if (detRs.size() != predictedComponents.size() || chi2s.size() != predictedComponents.size()) {
0083 edm::LogError("PosteriorWeightsCalculator") << "Problem in vector sizes";
0084 return std::vector<double>();
0085 }
0086
0087
0088
0089
0090
0091
0092 double sumWeights(0.);
0093 for (unsigned int i = 0; i < predictedComponents.size(); i++) {
0094 double priorWeight = predictedComponents[i].weight();
0095
0096 double chi2 = chi2s[i] - chi2Min;
0097
0098 double tempWeight(0.);
0099 if (detRs[i] > FLT_MIN) {
0100
0101
0102
0103
0104 tempWeight = priorWeight * std::sqrt(1. / detRs[i]) * std::exp(-0.5 * chi2);
0105 } else {
0106 LogDebug("GsfTrackFitters") << "PosteriorWeightsCalculator: detR < FLT_MIN !!";
0107 }
0108 weights.push_back(tempWeight);
0109 sumWeights += tempWeight;
0110 }
0111
0112 if (sumWeights < DBL_MIN) {
0113 LogDebug("GsfTrackFitters") << "PosteriorWeightsCalculator: sumWeight < DBL_MIN";
0114 edm::LogError("PosteriorWeightsCalculator") << "PosteriorWeightsCalculator: sumWeight < DBL_MIN";
0115 return std::vector<double>();
0116 }
0117
0118 if (weights.size() != predictedComponents.size()) {
0119 edm::LogError("PosteriorWeightsCalculator") << "Problem in vector sizes (2)";
0120 return std::vector<double>();
0121 }
0122 sumWeights = 1. / sumWeights;
0123 for (auto& w : weights)
0124 w *= sumWeights;
0125 return weights;
0126 }