TrackingTools/src/MuonBestMeasurementFinder.cc

0001 /** \class MuonBestMeasurementFinder
0002  *  Algorithmic class to get best measurement from a list of TM
0003  *  the chi2 cut for the MeasurementEstimator is huge since should not be used.
0004  *  The aim of this class is to return the "best" measurement according to the
0005  *  chi2, but without any cut. The decision whether to use or not the
0006  *  measurement is taken in the caller class.
0007  *  The evaluation is made (in hard-code way) with the granularity = 1. Where
0008  *  the granularity is the one defined in the MuonTrajectoyUpdatorClass.
0009  *
0010  *  \author R. Bellan - INFN Torino <riccardo.bellan@cern.ch>
0011  *  \author S. Lacaprara - INFN Legnaro <stefano.lacaprara@pd.infn.it>
0012  */
0013
0014 #include "RecoMuon/TrackingTools/interface/MuonBestMeasurementFinder.h"
0015
0016 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
0017 #include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHit.h"
0018
0019 #include "TrackingTools/PatternTools/interface/TrajectoryMeasurement.h"
0020 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h"
0021 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0022 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"
0023
0024 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0025
0026 using namespace std;
0027
0028 MuonBestMeasurementFinder::MuonBestMeasurementFinder() { theEstimator = new Chi2MeasurementEstimator(100000.); }
0029
0030 MuonBestMeasurementFinder::~MuonBestMeasurementFinder() { delete theEstimator; }
0031
0032 TrajectoryMeasurement* MuonBestMeasurementFinder::findBestMeasurement(std::vector<TrajectoryMeasurement>& measC,
0033                                                                       const Propagator* propagator) {
0034   const std::string metname = "Muon|RecoMuon|MuonBestMeasurementFinder";
0035
0036   TMContainer validMeasurements;
0037
0038   TrajectoryMeasurement* bestMeasurement = nullptr;
0039
0040   // consider only valid TM
0041   int NumValidMeas = 0;
0042   for (vector<TrajectoryMeasurement>::iterator measurement = measC.begin(); measurement != measC.end(); ++measurement) {
0043     if ((*measurement).recHit()->isValid()) {
0044       ++NumValidMeas;
0045       bestMeasurement = &(*measurement);
0046       validMeasurements.push_back(&(*measurement));
0047     }
0048   }
0049
0050   // If we have just one (or zero) valid meas, return it at once
0051   // (or return null measurement)
0052   if (NumValidMeas <= 1) {
0053     LogTrace(metname) << "MuonBestMeasurement: just " << NumValidMeas << " valid measurement ";
0054     return bestMeasurement;
0055   }
0056
0057   TMIterator measurement;
0058   double minChi2PerNDoF = 1.E6;
0059
0060   // if there are more than one valid measurement, then sort them.
0061   for (measurement = validMeasurements.begin(); measurement != validMeasurements.end(); measurement++) {
0062     TransientTrackingRecHit::ConstRecHitPointer muonRecHit = (*measurement)->recHit();
0063
0064     // FIXME!! FIXME !! FIXME !!
0065     pair<double, int> chi2Info = lookAtSubRecHits(*measurement, propagator);
0066
0067     double chi2PerNDoF = chi2Info.first / chi2Info.second;
0068     LogTrace(metname) << " The measurement has a chi2/npts " << chi2PerNDoF << " with dof = " << chi2Info.second
0069                       << " \n Till now the best chi2 is " << minChi2PerNDoF;
0070
0071     if (chi2PerNDoF && chi2PerNDoF < minChi2PerNDoF) {
0072       minChi2PerNDoF = chi2PerNDoF;
0073       bestMeasurement = *measurement;
0074     }
0075   }
0076   LogTrace(metname) << "The final best chi2 is " << minChi2PerNDoF << endl;
0077   return bestMeasurement;
0078 }
0079
0080 pair<double, int> MuonBestMeasurementFinder::lookAtSubRecHits(TrajectoryMeasurement* measurement,
0081                                                               const Propagator* propagator) {
0082   const std::string metname = "Muon|RecoMuon|MuonBestMeasurementFinder";
0083
0084   unsigned int npts = 0;
0085   // unused  double thisChi2 = 0.;
0086
0087   const TransientTrackingRecHit::ConstRecHitPointer& muonRecHit = measurement->recHit();
0088   const TrajectoryStateOnSurface& predState = measurement->predictedState();  // temporarily introduced by DT
0089   std::pair<bool, double> sing_chi2 =
0090       estimator()->estimate(predState, *(muonRecHit.get()));                  // temporarily introduced by DT
0091   npts = 1;                                                                   // temporarily introduced by DT
0092   std::pair<double, int> result = pair<double, int>(sing_chi2.second, npts);  // temporarily introduced by DT
0093
0094   //   // ask for the 2D-segments/2D-rechit                                                      // temporarily excluded by DT
0095   //   TransientTrackingRecHit::ConstRecHitContainer rhits_list = muonRecHit->transientHits();
0096
0097   //   LogTrace(metname)<<"Number of rechits in the measurement rechit: "<<rhits_list.size()<<endl;
0098
0099   //   // loop over them
0100   //   for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator rhit = rhits_list.begin();
0101   //        rhit!= rhits_list.end(); ++rhit)
0102   //     if ((*rhit)->isValid() ) {
0103   //       LogTrace(metname)<<"Rechit dimension: "<<(*rhit)->dimension()<<endl;
0104   //       npts+=(*rhit)->dimension();
0105
0106   //       TrajectoryStateOnSurface predState;
0107
0108   //       if (!( (*rhit)->geographicalId() == muonRecHit->geographicalId() ) ){
0109   //    predState = propagator->propagate(*measurement->predictedState().freeState(),
0110   //                      (*rhit)->det()->surface());
0111   //       }
0112   //       else predState = measurement->predictedState();
0113
0114   //       if ( predState.isValid() ) {
0115   //    std::pair<bool,double> sing_chi2 = estimator()->estimate( predState, *((*rhit).get()));
0116   //    thisChi2 += sing_chi2.second ;
0117   //    LogTrace(metname) << " single incremental chi2: " << sing_chi2.second;
0118   //       }
0119   //     }
0120
0121   //   pair<double,int> result = pair<double,int>(thisChi2,npts);
0122   return result;
0123 }