TrackingTools/src/MuonTrajectoryUpdator.cc

0001 /** \class MuonTrajectoryUpdator
0002  *  An updator for the Muon system
0003  *  This class update a trajectory with a muon chamber measurement.
0004  *  In spite of the name, it is NOT an updator, but has one.
0005  *  A muon RecHit is a segment (for DT and CSC) or a "hit" (RPC).
0006  *  This updator is suitable both for FW and BW filtering. The difference between the two fitter are two:
0007  *  the granularity of the updating (i.e.: segment position or 1D rechit position), which can be set via
0008  *  parameter set, and the propagation direction which is embeded in the propagator set in the c'tor.
0009  *
0010  *  \author R. Bellan - INFN Torino <riccardo.bellan@cern.ch>
0011  *  \author S. Lacaprara - INFN Legnaro
0012  *
0013  *  Modified by C. Calabria: GEM Implementation
0014  *  Modified by D. Nash: ME0 Implementation
0015  */
0016
0017 #include "RecoMuon/TrackingTools/interface/MuonTrajectoryUpdator.h"
0018 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
0019 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
0020
0021 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"
0022 #include "TrackingTools/PatternTools/interface/TrajectoryMeasurement.h"
0023 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0024 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0025 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0026 #include "TrackingTools/KalmanUpdators/interface/KFUpdator.h"
0027 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0028
0029 #include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHit.h"
0030 #include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHitBreaker.h"
0031
0032 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0033 #include "FWCore/Utilities/interface/Exception.h"
0034 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0035
0036 #include <algorithm>
0037
0038 using namespace edm;
0039 using namespace std;
0040
0041 /// Constructor from Propagator and Parameter set
0042 MuonTrajectoryUpdator::MuonTrajectoryUpdator(const edm::ParameterSet& par, NavigationDirection fitDirection)
0043     : theFitDirection(fitDirection) {
0044   // The max allowed chi2 to accept a rechit in the fit
0045   theMaxChi2 = par.getParameter<double>("MaxChi2");
0046   theEstimator = new Chi2MeasurementEstimator(theMaxChi2);
0047
0048   // The KF updator
0049   theUpdator = new KFUpdator();
0050
0051   // The granularity
0052   theGranularity = par.getParameter<int>("Granularity");
0053
0054   // Rescale the error of the first state?
0055   theRescaleErrorFlag = par.getParameter<bool>("RescaleError");
0056
0057   if (theRescaleErrorFlag)
0058     // The rescale factor
0059     theRescaleFactor = par.getParameter<double>("RescaleErrorFactor");
0060
0061   // Use the invalid hits?
0062   useInvalidHits = par.getParameter<bool>("UseInvalidHits");
0063
0064   // Flag needed for the rescaling
0065   theFirstTSOSFlag = true;
0066
0067   // Exlude RPC from the fit?
0068   theRPCExFlag = par.getParameter<bool>("ExcludeRPCFromFit");
0069 }
0070
0071 MuonTrajectoryUpdator::MuonTrajectoryUpdator(NavigationDirection fitDirection, double chi2, int granularity)
0072     : theMaxChi2(chi2), theGranularity(granularity), theFitDirection(fitDirection) {
0073   theEstimator = new Chi2MeasurementEstimator(theMaxChi2);
0074
0075   // The KF updator
0076   theUpdator = new KFUpdator();
0077 }
0078
0079 /// Destructor
0080 MuonTrajectoryUpdator::~MuonTrajectoryUpdator() {
0081   delete theEstimator;
0082   delete theUpdator;
0083 }
0084
0085 void MuonTrajectoryUpdator::makeFirstTime() { theFirstTSOSFlag = true; }
0086
0087 pair<bool, TrajectoryStateOnSurface> MuonTrajectoryUpdator::update(const TrajectoryMeasurement* measurement,
0088                                                                    Trajectory& trajectory,
0089                                                                    const Propagator* propagator) {
0090   const std::string metname = "Muon|RecoMuon|MuonTrajectoryUpdator";
0091
0092   MuonPatternRecoDumper muonDumper;
0093
0094   // Status of the updating
0095   bool updated = false;
0096
0097   if (!measurement)
0098     return pair<bool, TrajectoryStateOnSurface>(updated, TrajectoryStateOnSurface());
0099
0100   // measurement layer
0101   const DetLayer* detLayer = measurement->layer();
0102
0103   // these are the 4D segment for the CSC/DT and a point for the RPC
0104   TransientTrackingRecHit::ConstRecHitPointer muonRecHit = measurement->recHit();
0105
0106   // The KFUpdator takes TransientTrackingRecHits as arg.
0107   TransientTrackingRecHit::ConstRecHitContainer recHitsForFit =
0108       MuonTransientTrackingRecHitBreaker::breakInSubRecHits(muonRecHit, theGranularity);
0109
0110   // sort the container in agreement with the porpagation direction
0111   sort(recHitsForFit, detLayer);
0112
0113   TrajectoryStateOnSurface lastUpdatedTSOS = measurement->predictedState();
0114
0115   LogTrace(metname) << "Number of rechits for the fit: " << recHitsForFit.size() << endl;
0116
0117   TransientTrackingRecHit::ConstRecHitContainer::iterator recHit;
0118   for (recHit = recHitsForFit.begin(); recHit != recHitsForFit.end(); ++recHit) {
0119     if ((*recHit)->isValid()) {
0120       // propagate the TSOS onto the rechit plane
0121       TrajectoryStateOnSurface propagatedTSOS = propagateState(lastUpdatedTSOS, measurement, *recHit, propagator);
0122
0123       if (propagatedTSOS.isValid()) {
0124         pair<bool, double> thisChi2 = estimator()->estimate(propagatedTSOS, *((*recHit).get()));
0125
0126         LogTrace(metname) << "Estimation for Kalman Fit. Chi2: " << thisChi2.second;
0127
0128         // Is an RPC hit? Prepare the variable to possibly exluding it from the fit
0129         bool wantIncludeThisHit = true;
0130         if (theRPCExFlag && (*recHit)->geographicalId().det() == DetId::Muon &&
0131             (*recHit)->geographicalId().subdetId() == MuonSubdetId::RPC) {
0132           wantIncludeThisHit = false;
0133           LogTrace(metname)
0134               << "This is an RPC hit and the present configuration is such that it will be excluded from the fit";
0135         }
0136
0137         // The Chi2 cut was already applied in the estimator, which
0138         // returns 0 if the chi2 is bigger than the cut defined in its
0139         // constructor
0140         if (thisChi2.first) {
0141           updated = true;
0142           if (wantIncludeThisHit) {  // This split is a trick to have the RPC hits counted as updatable (in used chamber counting), while are not actually included in the fit when the proper obtion is activated.
0143
0144             LogTrace(metname) << endl
0145                               << "     Kalman Start"
0146                               << "\n"
0147                               << "\n";
0148             LogTrace(metname) << "  Meas. Position : " << (**recHit).globalPosition() << "\n"
0149                               << "  Pred. Position : " << propagatedTSOS.globalPosition()
0150                               << "  Pred Direction : " << propagatedTSOS.globalDirection() << endl;
0151
0152             if (theRescaleErrorFlag && theFirstTSOSFlag) {
0153               propagatedTSOS.rescaleError(theRescaleFactor);
0154               theFirstTSOSFlag = false;
0155             }
0156
0157             lastUpdatedTSOS = measurementUpdator()->update(propagatedTSOS, *((*recHit).get()));
0158
0159             LogTrace(metname) << "  Fit   Position : " << lastUpdatedTSOS.globalPosition()
0160                               << "  Fit  Direction : " << lastUpdatedTSOS.globalDirection() << "\n"
0161                               << "  Fit position radius : " << lastUpdatedTSOS.globalPosition().perp()
0162                               << "filter updated" << endl;
0163
0164             LogTrace(metname) << muonDumper.dumpTSOS(lastUpdatedTSOS);
0165
0166             LogTrace(metname) << "\n\n     Kalman End"
0167                               << "\n"
0168                               << "\n";
0169
0170             TrajectoryMeasurement&& updatedMeasurement =
0171                 updateMeasurement(propagatedTSOS, lastUpdatedTSOS, *recHit, thisChi2.second, detLayer, measurement);
0172             // FIXME: check!
0173             trajectory.push(std::move(updatedMeasurement), thisChi2.second);
0174           } else {
0175             LogTrace(metname) << "  Compatible RecHit with good chi2 but made with RPC when it was decided to not "
0176                                  "include it in the fit"
0177                               << "  --> trajectory NOT updated, invalid RecHit added." << endl;
0178
0179             MuonTransientTrackingRecHit::MuonRecHitPointer invalidRhPtr =
0180                 MuonTransientTrackingRecHit::specificBuild((*recHit)->det(), (*recHit)->hit());
0181             invalidRhPtr->invalidateHit();
0182             TrajectoryMeasurement invalidRhMeasurement(
0183                 propagatedTSOS, propagatedTSOS, invalidRhPtr, thisChi2.second, detLayer);
0184             trajectory.push(std::move(invalidRhMeasurement), thisChi2.second);
0185           }
0186         } else {
0187           if (useInvalidHits) {
0188             LogTrace(metname) << "  Compatible RecHit with too large chi2"
0189                               << "  --> trajectory NOT updated, invalid RecHit added." << endl;
0190
0191             MuonTransientTrackingRecHit::MuonRecHitPointer invalidRhPtr =
0192                 MuonTransientTrackingRecHit::specificBuild((*recHit)->det(), (*recHit)->hit());
0193             invalidRhPtr->invalidateHit();
0194             TrajectoryMeasurement invalidRhMeasurement(
0195                 propagatedTSOS, propagatedTSOS, invalidRhPtr, thisChi2.second, detLayer);
0196             trajectory.push(std::move(invalidRhMeasurement), thisChi2.second);
0197           }
0198         }
0199       }
0200     }
0201   }
0202   recHitsForFit.clear();
0203   return pair<bool, TrajectoryStateOnSurface>(updated, lastUpdatedTSOS);
0204 }
0205
0206 TrajectoryStateOnSurface MuonTrajectoryUpdator::propagateState(
0207     const TrajectoryStateOnSurface& state,
0208     const TrajectoryMeasurement* measurement,
0209     const TransientTrackingRecHit::ConstRecHitPointer& current,
0210     const Propagator* propagator) const {
0211   const TransientTrackingRecHit::ConstRecHitPointer& mother = measurement->recHit();
0212
0213   if (current->geographicalId() == mother->geographicalId())
0214     return measurement->predictedState();
0215
0216   const TrajectoryStateOnSurface tsos = propagator->propagate(state, current->det()->surface());
0217   return tsos;
0218 }
0219
0220 // FIXME: would I a different threatment for the two prop dirrections??
0221 TrajectoryMeasurement MuonTrajectoryUpdator::updateMeasurement(const TrajectoryStateOnSurface& propagatedTSOS,
0222                                                                const TrajectoryStateOnSurface& lastUpdatedTSOS,
0223                                                                const TransientTrackingRecHit::ConstRecHitPointer& recHit,
0224                                                                const double& chi2,
0225                                                                const DetLayer* detLayer,
0226                                                                const TrajectoryMeasurement* initialMeasurement) {
0227   return TrajectoryMeasurement(propagatedTSOS, lastUpdatedTSOS, recHit, chi2, detLayer);
0228
0229   //   // FIXME: put a better check! One could fit in first out-in and then in - out
0230   //   if(propagator()->propagationDirection() == alongMomentum)
0231   //     return TrajectoryMeasurement(propagatedTSOS, lastUpdatedTSOS,
0232   //                 recHit,thisChi2.second,detLayer);
0233
0234   //   // FIXME: Check this carefully!!
0235   //   else if(propagator()->propagationDirection() == oppositeToMomentum)
0236   //     return TrajectoryMeasurement(initialMeasurement->forwardPredictedState(),
0237   //                 propagatedTSOS, lastUpdatedTSOS,
0238   //                 recHit,thisChi2.second,detLayer);
0239   //   else{
0240   //     LogError("MuonTrajectoryUpdator::updateMeasurement") <<"Wrong propagation direction!!";
0241   //   }
0242 }
0243
0244 void MuonTrajectoryUpdator::sort(TransientTrackingRecHit::ConstRecHitContainer& recHitsForFit,
0245                                  const DetLayer* detLayer) {
0246   if (detLayer->subDetector() == GeomDetEnumerators::DT) {
0247     if (fitDirection() == insideOut)
0248       stable_sort(recHitsForFit.begin(), recHitsForFit.end(), RadiusComparatorInOut());
0249     else if (fitDirection() == outsideIn)
0250       stable_sort(recHitsForFit.begin(), recHitsForFit.end(), RadiusComparatorOutIn());
0251     else
0252       LogError("Muon|RecoMuon|MuonTrajectoryUpdator") << "MuonTrajectoryUpdator::sort: Wrong propagation direction!!";
0253   }
0254
0255   else if (detLayer->subDetector() == GeomDetEnumerators::CSC) {
0256     if (fitDirection() == insideOut)
0257       stable_sort(recHitsForFit.begin(), recHitsForFit.end(), ZedComparatorInOut());
0258     else if (fitDirection() == outsideIn)
0259       stable_sort(recHitsForFit.begin(), recHitsForFit.end(), ZedComparatorOutIn());
0260     else
0261       LogError("Muon|RecoMuon|MuonTrajectoryUpdator") << "MuonTrajectoryUpdator::sort: Wrong propagation direction!!";
0262   }
0263
0264   else if (detLayer->subDetector() == GeomDetEnumerators::GEM) {
0265     if (fitDirection() == insideOut)
0266       stable_sort(recHitsForFit.begin(), recHitsForFit.end(), ZedComparatorInOut());
0267     else if (fitDirection() == outsideIn)
0268       stable_sort(recHitsForFit.begin(), recHitsForFit.end(), ZedComparatorOutIn());
0269     else
0270       LogError("Muon|RecoMuon|MuonTrajectoryUpdator") << "MuonTrajectoryUpdator::sort: Wrong propagation direction!!";
0271   } else if (detLayer->subDetector() == GeomDetEnumerators::ME0) {
0272     if (fitDirection() == insideOut)
0273       stable_sort(recHitsForFit.begin(), recHitsForFit.end(), ZedComparatorInOut());
0274     else if (fitDirection() == outsideIn)
0275       stable_sort(recHitsForFit.begin(), recHitsForFit.end(), ZedComparatorOutIn());
0276     else
0277       LogError("Muon|RecoMuon|MuonTrajectoryUpdator") << "MuonTrajectoryUpdator::sort: Wrong propagation direction!!";
0278   }
0279 }