GsfTracking/src/PosteriorWeightsCalculator.cc

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   // calculate chi2 and determinant / component and find
0052   //   minimum / maximum of chi2
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   // calculate weights (extracting a common factor
0089   //   exp(-0.5*chi2Min) to avoid numerical problems
0090   //   during exponentation
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       // Calculation of (non-normalised) weight. Common factors exp(-chi2Norm/2.) and
0102       // 1./sqrt(2*pi*recHit.dimension()) have been omitted
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 }