KalmanUpdators/src/MRHChi2MeasurementEstimator.cc

0001 #include "DataFormats/TrackerRecHit2D/interface/SiTrackerMultiRecHit.h"
0002 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0003 #include "TrackingTools/KalmanUpdators/interface/MRHChi2MeasurementEstimator.h"
0004 #include "TrackingTools/PatternTools/interface/MeasurementExtractor.h"
0005 #include "DataFormats/GeometrySurface/interface/Plane.h"
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 #include "DataFormats/Math/interface/invertPosDefMatrix.h"
0008 #include "DataFormats/Math/interface/ProjectMatrix.h"
0009 #include "DataFormats/TrackingRecHit/interface/KfComponentsHolder.h"
0010
0011 std::pair<bool, double> MRHChi2MeasurementEstimator::estimate(const TrajectoryStateOnSurface& tsos,
0012                                                               const TrackingRecHit& aRecHit) const {
0013   switch (aRecHit.dimension()) {
0014     case 1:
0015       return estimate<1>(tsos, aRecHit);
0016     case 2:
0017       return estimate<2>(tsos, aRecHit);
0018     //avoid the not-(1D or 2D)  hit due to the final sum
0019     //supposing all the hits inside of a MRH have the same dimension
0020     case 3:
0021     case 4:
0022     case 5: {
0023       LogDebug("MRHChi2MeasurementEstimator") << "WARNING:The hit is not 1D either 2D:"
0024                                               << " does not count in the MRH Chi2 estimation.";
0025       double est = 0.0;
0026       return HitReturnType(false, est);
0027     }
0028   }
0029   throw cms::Exception("Rec hit of invalid dimension (not 1,2,3,4,5)")
0030       << "The value was " << aRecHit.dimension() << ", type is " << typeid(aRecHit).name() << "\n";
0031 }
0032
0033 //---------------------------------------------------------------------------------------------------------------
0034 template <unsigned int N>
0035 std::pair<bool, double> MRHChi2MeasurementEstimator::estimate(const TrajectoryStateOnSurface& tsos,
0036                                                               const TrackingRecHit& aRecHit) const {
0037   LogDebug("MRHChi2MeasurementEstimator") << "Calling MRHChi2MeasurementEstimator";
0038   SiTrackerMultiRecHit const& mHit = dynamic_cast<SiTrackerMultiRecHit const&>(aRecHit);
0039   double est = 0;
0040
0041   double annealing = mHit.getAnnealingFactor();
0042   LogDebug("MRHChi2MeasurementEstimator") << "Current annealing factor is " << annealing;
0043
0044   std::vector<const TrackingRecHit*> components = mHit.recHits();
0045   LogDebug("MRHChi2MeasurementEstimator") << "this hit has " << components.size() << " components";
0046
0047   int iComp = 0;
0048   for (std::vector<const TrackingRecHit*>::const_iterator iter = components.begin(); iter != components.end();
0049        iter++, iComp++) {
0050     // define variables that will be used to setup the KfComponentsHolder
0051     ProjectMatrix<double, 5, N> pf;
0052     typename AlgebraicROOTObject<N>::Vector r, rMeas;
0053     typename AlgebraicROOTObject<N, N>::SymMatrix V, VMeas;
0054     AlgebraicVector5 x = tsos.localParameters().vector();
0055     const AlgebraicSymMatrix55& C = (tsos.localError().matrix());
0056
0057     // setup the holder with the correct dimensions and get the values
0058     KfComponentsHolder holder;
0059     holder.template setup<N>(&r, &V, &pf, &rMeas, &VMeas, x, C);
0060     (**iter).getKfComponents(holder);
0061
0062     r -= rMeas;
0063     V = V * annealing + VMeas;
0064     bool ierr = invertPosDefMatrix(V);
0065     if (!ierr) {
0066       edm::LogError("SiTrackerMultiRecHitUpdator")
0067           << "SiTrackerMultiRecHitUpdator::ComputeParameters2dim: W not valid!" << std::endl;
0068     }
0069
0070     LogDebug("MRHChi2MeasurementEstimator") << "Hit with weight " << mHit.weight(iComp);
0071     est += ROOT::Math::Similarity(r, V) * mHit.weight(iComp);
0072   }
0073
0074   return returnIt(est);
0075 }