StandAloneTrackFinder/src/StandAloneMuonFilter.cc

0001 /** \class StandAloneMuonFilter
0002  *  The inward-outward fitter (starts from seed state).
0003  *
0004  *  \author R. Bellan - INFN Torino <riccardo.bellan@cern.ch>
0005  *          D. Trocino - INFN Torino <daniele.trocino@to.infn.it>
0006  *
0007  *  Modified by C. Calabria
0008  *  Modified by D. Nash
0009  */
0010 #include "RecoMuon/StandAloneTrackFinder/interface/StandAloneMuonFilter.h"
0011
0012 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0013
0014 // FIXME: remove this
0015 #include "FWCore/Framework/interface/Event.h"
0016
0017 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
0018 #include "RecoMuon/TrackingTools/interface/MuonTrajectoryUpdator.h"
0019 #include "RecoMuon/MeasurementDet/interface/MuonDetLayerMeasurements.h"
0020 #include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
0021 #include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHit.h"
0022 #include "RecoMuon/Navigation/interface/DirectMuonNavigation.h"
0023
0024 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0025
0026 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"
0027
0028 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0029
0030 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0031 #include "TrackingTools/PatternTools/interface/TrajectoryMeasurement.h"
0032 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0033
0034 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0035 #include "FWCore/Utilities/interface/Exception.h"
0036
0037 #include <vector>
0038
0039 using namespace edm;
0040 using namespace std;
0041
0042 StandAloneMuonFilter::StandAloneMuonFilter(const ParameterSet& par,
0043                                            const MuonServiceProxy* service,
0044                                            edm::ConsumesCollector& iC)
0045     : theService(service), theOverlappingChambersFlag(true) {
0046   // Fit direction
0047   string fitDirectionName = par.getParameter<string>("FitDirection");
0048
0049   if (fitDirectionName == "insideOut")
0050     theFitDirection = insideOut;
0051   else if (fitDirectionName == "outsideIn")
0052     theFitDirection = outsideIn;
0053   else
0054     throw cms::Exception("StandAloneMuonFilter constructor")
0055         << "Wrong fit direction chosen in StandAloneMuonFilter::StandAloneMuonFilter ParameterSet"
0056         << "\n"
0057         << "Possible choices are:"
0058         << "\n"
0059         << "FitDirection = insideOut or FitDirection = outsideIn";
0060
0061   // The max allowed chi2 to accept a rechit in the fit
0062   theMaxChi2 = par.getParameter<double>("MaxChi2");
0063
0064   // The errors of the trajectory state are multiplied by nSigma
0065   // to define acceptance of BoundPlane and maximalLocalDisplacement
0066   theNSigma = par.getParameter<double>("NumberOfSigma");  // default = 3.
0067
0068   // The navigation type:
0069   // "Direct","Standard"
0070   theNavigationType = par.getParameter<string>("NavigationType");
0071
0072   // The estimator: makes the decision wheter a measure is good or not
0073   // it isn't used by the updator which does the real fit. In fact, in principle,
0074   // a looser request onto the measure set can be requested
0075   // (w.r.t. the request on the accept/reject measure in the fit)
0076   theEstimator = new Chi2MeasurementEstimator(theMaxChi2, theNSigma);
0077
0078   thePropagatorName = par.getParameter<string>("Propagator");
0079
0080   theBestMeasurementFinder = new MuonBestMeasurementFinder();
0081
0082   // Muon trajectory updator parameters
0083   ParameterSet muonUpdatorPSet = par.getParameter<ParameterSet>("MuonTrajectoryUpdatorParameters");
0084
0085   // the updator needs the fit direction
0086   theMuonUpdator = new MuonTrajectoryUpdator(muonUpdatorPSet, fitDirection());
0087
0088   // Measurement Extractor: enable the measure for each muon sub detector
0089   bool enableDTMeasurement = par.getParameter<bool>("EnableDTMeasurement");
0090   bool enableCSCMeasurement = par.getParameter<bool>("EnableCSCMeasurement");
0091   bool enableRPCMeasurement = par.getParameter<bool>("EnableRPCMeasurement");
0092   bool enableGEMMeasurement = par.getParameter<bool>("EnableGEMMeasurement");
0093   bool enableME0Measurement = par.getParameter<bool>("EnableME0Measurement");
0094
0095   theMeasurementExtractor = new MuonDetLayerMeasurements(par.getParameter<InputTag>("DTRecSegmentLabel"),
0096                                                          par.getParameter<InputTag>("CSCRecSegmentLabel"),
0097                                                          par.getParameter<InputTag>("RPCRecSegmentLabel"),
0098                                                          par.getParameter<InputTag>("GEMRecSegmentLabel"),
0099                                                          par.getParameter<InputTag>("ME0RecSegmentLabel"),
0100                                                          iC,
0101                                                          enableDTMeasurement,
0102                                                          enableCSCMeasurement,
0103                                                          enableRPCMeasurement,
0104                                                          enableGEMMeasurement,
0105                                                          enableME0Measurement);
0106
0107   theRPCLoneliness = (!(enableDTMeasurement && enableCSCMeasurement)) ? enableRPCMeasurement : false;
0108 }
0109
0110 StandAloneMuonFilter::~StandAloneMuonFilter() {
0111   LogTrace("Muon|RecoMuon|StandAloneMuonFilter") << "StandAloneMuonFilter destructor called" << endl;
0112
0113   delete theEstimator;
0114   delete theMuonUpdator;
0115   delete theMeasurementExtractor;
0116   delete theBestMeasurementFinder;
0117 }
0118
0119 const Propagator* StandAloneMuonFilter::propagator() const { return &*theService->propagator(thePropagatorName); }
0120
0121 /// Return the propagation direction
0122 PropagationDirection StandAloneMuonFilter::propagationDirection() const {
0123   if (fitDirection() == 0)
0124     return alongMomentum;
0125   else if (fitDirection() == 1)
0126     return oppositeToMomentum;
0127   else
0128     return anyDirection;
0129 }
0130
0131 void StandAloneMuonFilter::reset() {
0132   totalChambers = dtChambers = cscChambers = rpcChambers = gemChambers = me0Chambers = 0;
0133   totalCompatibleChambers = dtCompatibleChambers = cscCompatibleChambers = rpcCompatibleChambers =
0134       gemCompatibleChambers = me0CompatibleChambers = 0;
0135
0136   theLastCompatibleTSOS = theLastUpdatedTSOS = theLastButOneUpdatedTSOS = TrajectoryStateOnSurface();
0137
0138   theMuonUpdator->makeFirstTime();
0139
0140   theDetLayers.clear();
0141 }
0142
0143 void StandAloneMuonFilter::setEvent(const Event& event) { theMeasurementExtractor->setEvent(event); }
0144
0145 void StandAloneMuonFilter::incrementChamberCounters(const DetLayer* layer) {
0146   if (layer->subDetector() == GeomDetEnumerators::DT)
0147     dtChambers++;
0148   else if (layer->subDetector() == GeomDetEnumerators::CSC)
0149     cscChambers++;
0150   else if (layer->subDetector() == GeomDetEnumerators::RPCBarrel ||
0151            layer->subDetector() == GeomDetEnumerators::RPCEndcap)
0152     rpcChambers++;
0153   else if (layer->subDetector() == GeomDetEnumerators::GEM)
0154     gemChambers++;
0155   else if (layer->subDetector() == GeomDetEnumerators::ME0)
0156     me0Chambers++;
0157   else
0158     LogError("Muon|RecoMuon|StandAloneMuonFilter") << "Unrecognized module type in incrementChamberCounters";
0159   // FIXME:
0160   //   << layer->module() << " " <<layer->Part() << endl;
0161
0162   totalChambers++;
0163 }
0164
0165 void StandAloneMuonFilter::incrementCompatibleChamberCounters(const DetLayer* layer) {
0166   if (layer->subDetector() == GeomDetEnumerators::DT)
0167     dtCompatibleChambers++;
0168   else if (layer->subDetector() == GeomDetEnumerators::CSC)
0169     cscCompatibleChambers++;
0170   else if (layer->subDetector() == GeomDetEnumerators::RPCBarrel ||
0171            layer->subDetector() == GeomDetEnumerators::RPCEndcap)
0172     rpcCompatibleChambers++;
0173   else if (layer->subDetector() == GeomDetEnumerators::GEM)
0174     gemCompatibleChambers++;
0175   else if (layer->subDetector() == GeomDetEnumerators::ME0)
0176     me0CompatibleChambers++;
0177   else
0178     LogError("Muon|RecoMuon|StandAloneMuonFilter") << "Unrecognized module type in incrementCompatibleChamberCounters";
0179
0180   totalCompatibleChambers++;
0181 }
0182
0183 vector<const DetLayer*> StandAloneMuonFilter::compatibleLayers(const DetLayer* initialLayer,
0184                                                                const FreeTrajectoryState& fts,
0185                                                                PropagationDirection propDir) {
0186   vector<const DetLayer*> detLayers;
0187
0188   if (theNavigationType == "Standard") {
0189     // ask for compatible layers
0190     detLayers = theService->muonNavigationSchool()->compatibleLayers(*initialLayer, fts, propDir);
0191     // I have to fit by hand the first layer until the seedTSOS is defined on the first rechit layer
0192     // In fact the first layer is not returned by initialLayer->compatibleLayers.
0193     detLayers.insert(detLayers.begin(), initialLayer);
0194   } else if (theNavigationType == "Direct") {
0195     DirectMuonNavigation navigation(theService->detLayerGeometry());
0196     detLayers = navigation.compatibleLayers(fts, propDir);
0197   } else
0198     edm::LogError("Muon|RecoMuon|StandAloneMuonFilter") << "No Properly Navigation Selected!!" << endl;
0199
0200   return detLayers;
0201 }
0202
0203 void StandAloneMuonFilter::refit(const TrajectoryStateOnSurface& initialTSOS,
0204                                  const DetLayer* initialLayer,
0205                                  Trajectory& trajectory) {
0206   const std::string metname = "Muon|RecoMuon|StandAloneMuonFilter";
0207
0208   // reset the refitter each seed refinement
0209   reset();
0210
0211   MuonPatternRecoDumper debug;
0212
0213   LogTrace(metname) << "Starting the refit" << endl;
0214
0215   // this is the most outward TSOS (updated or predicted) onto a DetLayer
0216   TrajectoryStateOnSurface lastTSOS = theLastCompatibleTSOS = theLastUpdatedTSOS = theLastButOneUpdatedTSOS =
0217       initialTSOS;
0218
0219   double eta0 = initialTSOS.freeTrajectoryState()->momentum().eta();
0220   vector<const DetLayer*> detLayers =
0221       compatibleLayers(initialLayer, *initialTSOS.freeTrajectoryState(), propagationDirection());
0222
0223   LogTrace(metname) << "compatible layers found: " << detLayers.size() << endl;
0224
0225   vector<const DetLayer*>::const_iterator layer;
0226
0227   // the layers are ordered in agreement with the fit/propagation direction
0228   for (layer = detLayers.begin(); layer != detLayers.end(); ++layer) {
0229     //    bool firstTime = true;
0230
0231     LogTrace(metname) << debug.dumpLayer(*layer);
0232
0233     LogTrace(metname) << "search Trajectory Measurement from: " << lastTSOS.globalPosition();
0234
0235     // pick the best measurement from each group
0236     std::vector<TrajectoryMeasurement> bestMeasurements{};
0237
0238     if (lastTSOS.isValid())
0239       bestMeasurements = findBestMeasurements(*layer, lastTSOS);
0240     else
0241       edm::LogInfo(metname) << "Invalid last TSOS, will not find best measurements ";
0242
0243     // RB: Different ways can be choosen if no bestMeasurement is available:
0244     // 1- check on lastTSOS-initialTSOS eta difference
0245     // 2- check on lastTSOS-lastButOneUpdatedTSOS eta difference
0246     // After this choice:
0247     // A- extract the measurements compatible with the initialTSOS (seed)
0248     // B- extract the measurements compatible with the lastButOneUpdatedTSOS
0249     // In ORCA the choice was 1A. Here I will try 1B and if it fail I'll try 1A
0250     // another possibility could be 2B and then 1A.
0251
0252     // if no measurement found and the current TSOS has an eta very different
0253     // wrt the initial one (i.e. seed), then try to find the measurements
0254     // according to the lastButOne FTS. (1B)
0255     double lastdEta = fabs(lastTSOS.freeTrajectoryState()->momentum().eta() - eta0);
0256     if (bestMeasurements.empty() && lastdEta > 0.1) {
0257       LogTrace(metname) << "No measurement and big eta variation wrt seed" << endl
0258                         << "trying with lastButOneUpdatedTSOS";
0259
0260       if (theLastButOneUpdatedTSOS.isValid())
0261         bestMeasurements = findBestMeasurements(*layer, theLastButOneUpdatedTSOS);
0262       else
0263         edm::LogInfo(metname) << "Invalid last but one updated TSOS, will not find best measurements ";
0264     }
0265
0266     //if no measurement found and the current FTS has an eta very different
0267     //wrt the initial one (i.e. seed), then try to find the measurements
0268     //according to the initial FTS. (1A)
0269     if (bestMeasurements.empty() && lastdEta > 0.1) {
0270       LogTrace(metname) << "No measurement and big eta variation wrt seed" << endl << "tryng with seed TSOS";
0271
0272       if (initialTSOS.isValid())
0273         bestMeasurements = findBestMeasurements(*layer, initialTSOS);
0274       else
0275         edm::LogInfo(metname) << "Invalid initial TSOS, will not find best measurements ";
0276     }
0277
0278     // FIXME: uncomment this line!!
0279     // if(!bestMeasurement && firstTime) break;
0280
0281     if (!bestMeasurements.empty()) {
0282       incrementCompatibleChamberCounters(*layer);
0283       bool added = false;
0284       for (std::vector<TrajectoryMeasurement>::const_iterator tmItr = bestMeasurements.begin();
0285            tmItr != bestMeasurements.end();
0286            ++tmItr) {
0287         added |= update(*layer, &(*tmItr), trajectory);
0288         lastTSOS = theLastUpdatedTSOS;
0289       }
0290       if (added) {
0291         incrementChamberCounters(*layer);
0292         theDetLayers.push_back(*layer);
0293       }
0294     }
0295     // SL in case no valid mesurement is found, still I want to use the predicted
0296     // state for the following measurement serches. I take the first in the
0297     // container. FIXME!!! I want to carefully check this!!!!!
0298     else {
0299       LogTrace(metname) << "No best measurement found" << endl;
0300       //       if (!theMeasurementCache.empty()){
0301       //    LogTrace(metname)<<"but the #of measurement is "<<theMeasurementCache.size()<<endl;
0302       //         lastTSOS = theMeasurementCache.front().predictedState();
0303       //       }
0304     }
0305   }  // loop over layers
0306 }
0307
0308 std::vector<TrajectoryMeasurement> StandAloneMuonFilter::findBestMeasurements(const DetLayer* layer,
0309                                                                               const TrajectoryStateOnSurface& tsos) {
0310   const std::string metname = "Muon|RecoMuon|StandAloneMuonFilter";
0311
0312   std::vector<TrajectoryMeasurement> result;
0313   std::vector<TrajectoryMeasurement> measurements;
0314
0315   if (theOverlappingChambersFlag && layer->hasGroups()) {
0316     std::vector<TrajectoryMeasurementGroup> measurementGroups =
0317         theMeasurementExtractor->groupedMeasurements(layer, tsos, *propagator(), *estimator());
0318
0319     if (theFitDirection == outsideIn) {
0320       LogTrace(metname) << "Reversing the order of groupedMeasurements as the direction of the fit is outside-in";
0321       reverse(measurementGroups.begin(), measurementGroups.end());
0322       // this should be fixed either in RecoMuon/MeasurementDet/MuonDetLayerMeasurements or
0323       // RecoMuon/DetLayers/MuRingForwardDoubleLayer
0324     }
0325
0326     for (std::vector<TrajectoryMeasurementGroup>::const_iterator tmGroupItr = measurementGroups.begin();
0327          tmGroupItr != measurementGroups.end();
0328          ++tmGroupItr) {
0329       measurements = tmGroupItr->measurements();
0330       LogTrace(metname) << "Number of Trajectory Measurement: " << measurements.size();
0331
0332       const TrajectoryMeasurement* bestMeasurement =
0333           bestMeasurementFinder()->findBestMeasurement(measurements, propagator());
0334
0335       if (bestMeasurement)
0336         result.push_back(*bestMeasurement);
0337     }
0338   } else {
0339     measurements = theMeasurementExtractor->measurements(layer, tsos, *propagator(), *estimator());
0340     LogTrace(metname) << "Number of Trajectory Measurement: " << measurements.size();
0341     const TrajectoryMeasurement* bestMeasurement =
0342         bestMeasurementFinder()->findBestMeasurement(measurements, propagator());
0343     if (bestMeasurement)
0344       result.push_back(*bestMeasurement);
0345   }
0346   return result;
0347 }
0348
0349 bool StandAloneMuonFilter::update(const DetLayer* layer, const TrajectoryMeasurement* meas, Trajectory& trajectory) {
0350   const std::string metname = "Muon|RecoMuon|StandAloneMuonFilter";
0351   MuonPatternRecoDumper debug;
0352
0353   LogTrace(metname) << "best measurement found"
0354                     << "\n"
0355                     << "updating the trajectory..." << endl;
0356   pair<bool, TrajectoryStateOnSurface> result = updator()->update(meas, trajectory, propagator());
0357   LogTrace(metname) << "trajectory updated: " << result.first << endl;
0358   LogTrace(metname) << debug.dumpTSOS(result.second);
0359
0360   if (result.first) {
0361     theLastButOneUpdatedTSOS = theLastUpdatedTSOS;
0362     theLastUpdatedTSOS = result.second;
0363   }
0364
0365   if (result.second.isValid())
0366     theLastCompatibleTSOS = result.second;
0367
0368   return result.first;
0369 }
0370
0371 void StandAloneMuonFilter::createDefaultTrajectory(const Trajectory& oldTraj, Trajectory& defTraj) {
0372   Trajectory::DataContainer const& oldMeas = oldTraj.measurements();
0373   defTraj.reserve(oldMeas.size());
0374
0375   for (Trajectory::DataContainer::const_iterator itm = oldMeas.begin(); itm != oldMeas.end(); itm++) {
0376     if (!(*itm).recHit()->isValid())
0377       defTraj.push(*itm, (*itm).estimate());
0378     else {
0379       MuonTransientTrackingRecHit::MuonRecHitPointer invRhPtr =
0380           MuonTransientTrackingRecHit::specificBuild((*itm).recHit()->det(), (*itm).recHit()->hit());
0381       invRhPtr->invalidateHit();
0382       TrajectoryMeasurement invRhMeas(
0383           (*itm).forwardPredictedState(), (*itm).updatedState(), invRhPtr, (*itm).estimate(), (*itm).layer());
0384       defTraj.push(std::move(invRhMeas), (*itm).estimate());
0385     }
0386
0387   }  // end for
0388 }