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             if (!lastUpdatedTSOS.isValid()) {
0160               edm::LogInfo(metname) << "Invalid last TSOS, will skip RecHit ";
0161               lastUpdatedTSOS = propagatedTSOS;  // Revert update
0162               continue;
0163             }
0164
0165             LogTrace(metname) << "  Fit   Position : " << lastUpdatedTSOS.globalPosition()
0166                               << "  Fit  Direction : " << lastUpdatedTSOS.globalDirection() << "\n"
0167                               << "  Fit position radius : " << lastUpdatedTSOS.globalPosition().perp()
0168                               << "filter updated" << endl;
0169
0170             LogTrace(metname) << muonDumper.dumpTSOS(lastUpdatedTSOS);
0171
0172             LogTrace(metname) << "\n\n     Kalman End"
0173                               << "\n"
0174                               << "\n";
0175
0176             TrajectoryMeasurement&& updatedMeasurement =
0177                 updateMeasurement(propagatedTSOS, lastUpdatedTSOS, *recHit, thisChi2.second, detLayer, measurement);
0178             // FIXME: check!
0179             trajectory.push(std::move(updatedMeasurement), thisChi2.second);
0180           } else {
0181             LogTrace(metname) << "  Compatible RecHit with good chi2 but made with RPC when it was decided to not "
0182                                  "include it in the fit"
0183                               << "  --> trajectory NOT updated, invalid RecHit added." << endl;
0184
0185             MuonTransientTrackingRecHit::MuonRecHitPointer invalidRhPtr =
0186                 MuonTransientTrackingRecHit::specificBuild((*recHit)->det(), (*recHit)->hit());
0187             invalidRhPtr->invalidateHit();
0188             TrajectoryMeasurement invalidRhMeasurement(
0189                 propagatedTSOS, propagatedTSOS, invalidRhPtr, thisChi2.second, detLayer);
0190             trajectory.push(std::move(invalidRhMeasurement), thisChi2.second);
0191           }
0192         } else {
0193           if (useInvalidHits) {
0194             LogTrace(metname) << "  Compatible RecHit with too large chi2"
0195                               << "  --> trajectory NOT updated, invalid RecHit added." << endl;
0196
0197             MuonTransientTrackingRecHit::MuonRecHitPointer invalidRhPtr =
0198                 MuonTransientTrackingRecHit::specificBuild((*recHit)->det(), (*recHit)->hit());
0199             invalidRhPtr->invalidateHit();
0200             TrajectoryMeasurement invalidRhMeasurement(
0201                 propagatedTSOS, propagatedTSOS, invalidRhPtr, thisChi2.second, detLayer);
0202             trajectory.push(std::move(invalidRhMeasurement), thisChi2.second);
0203           }
0204         }
0205       }
0206     }
0207   }
0208   recHitsForFit.clear();
0209   return pair<bool, TrajectoryStateOnSurface>(updated, lastUpdatedTSOS);
0210 }
0211
0212 TrajectoryStateOnSurface MuonTrajectoryUpdator::propagateState(
0213     const TrajectoryStateOnSurface& state,
0214     const TrajectoryMeasurement* measurement,
0215     const TransientTrackingRecHit::ConstRecHitPointer& current,
0216     const Propagator* propagator) const {
0217   const TransientTrackingRecHit::ConstRecHitPointer& mother = measurement->recHit();
0218
0219   if (current->geographicalId() == mother->geographicalId())
0220     return measurement->predictedState();
0221
0222   const TrajectoryStateOnSurface tsos = propagator->propagate(state, current->det()->surface());
0223   return tsos;
0224 }
0225
0226 // FIXME: would I a different threatment for the two prop dirrections??
0227 TrajectoryMeasurement MuonTrajectoryUpdator::updateMeasurement(const TrajectoryStateOnSurface& propagatedTSOS,
0228                                                                const TrajectoryStateOnSurface& lastUpdatedTSOS,
0229                                                                const TransientTrackingRecHit::ConstRecHitPointer& recHit,
0230                                                                const double& chi2,
0231                                                                const DetLayer* detLayer,
0232                                                                const TrajectoryMeasurement* initialMeasurement) {
0233   return TrajectoryMeasurement(propagatedTSOS, lastUpdatedTSOS, recHit, chi2, detLayer);
0234
0235   //   // FIXME: put a better check! One could fit in first out-in and then in - out
0236   //   if(propagator()->propagationDirection() == alongMomentum)
0237   //     return TrajectoryMeasurement(propagatedTSOS, lastUpdatedTSOS,
0238   //                 recHit,thisChi2.second,detLayer);
0239
0240   //   // FIXME: Check this carefully!!
0241   //   else if(propagator()->propagationDirection() == oppositeToMomentum)
0242   //     return TrajectoryMeasurement(initialMeasurement->forwardPredictedState(),
0243   //                 propagatedTSOS, lastUpdatedTSOS,
0244   //                 recHit,thisChi2.second,detLayer);
0245   //   else{
0246   //     LogError("MuonTrajectoryUpdator::updateMeasurement") <<"Wrong propagation direction!!";
0247   //   }
0248 }
0249
0250 void MuonTrajectoryUpdator::sort(TransientTrackingRecHit::ConstRecHitContainer& recHitsForFit,
0251                                  const DetLayer* detLayer) {
0252   if (detLayer->subDetector() == GeomDetEnumerators::DT) {
0253     if (fitDirection() == insideOut)
0254       stable_sort(recHitsForFit.begin(), recHitsForFit.end(), RadiusComparatorInOut());
0255     else if (fitDirection() == outsideIn)
0256       stable_sort(recHitsForFit.begin(), recHitsForFit.end(), RadiusComparatorOutIn());
0257     else
0258       LogError("Muon|RecoMuon|MuonTrajectoryUpdator") << "MuonTrajectoryUpdator::sort: Wrong propagation direction!!";
0259   }
0260
0261   else if (detLayer->subDetector() == GeomDetEnumerators::CSC) {
0262     if (fitDirection() == insideOut)
0263       stable_sort(recHitsForFit.begin(), recHitsForFit.end(), ZedComparatorInOut());
0264     else if (fitDirection() == outsideIn)
0265       stable_sort(recHitsForFit.begin(), recHitsForFit.end(), ZedComparatorOutIn());
0266     else
0267       LogError("Muon|RecoMuon|MuonTrajectoryUpdator") << "MuonTrajectoryUpdator::sort: Wrong propagation direction!!";
0268   }
0269
0270   else if (detLayer->subDetector() == GeomDetEnumerators::GEM) {
0271     if (fitDirection() == insideOut)
0272       stable_sort(recHitsForFit.begin(), recHitsForFit.end(), ZedComparatorInOut());
0273     else if (fitDirection() == outsideIn)
0274       stable_sort(recHitsForFit.begin(), recHitsForFit.end(), ZedComparatorOutIn());
0275     else
0276       LogError("Muon|RecoMuon|MuonTrajectoryUpdator") << "MuonTrajectoryUpdator::sort: Wrong propagation direction!!";
0277   } else if (detLayer->subDetector() == GeomDetEnumerators::ME0) {
0278     if (fitDirection() == insideOut)
0279       stable_sort(recHitsForFit.begin(), recHitsForFit.end(), ZedComparatorInOut());
0280     else if (fitDirection() == outsideIn)
0281       stable_sort(recHitsForFit.begin(), recHitsForFit.end(), ZedComparatorOutIn());
0282     else
0283       LogError("Muon|RecoMuon|MuonTrajectoryUpdator") << "MuonTrajectoryUpdator::sort: Wrong propagation direction!!";
0284   }
0285 }