Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:26:13

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 = findBestMeasurements(*layer, lastTSOS);
0237 
0238     // RB: Different ways can be choosen if no bestMeasurement is available:
0239     // 1- check on lastTSOS-initialTSOS eta difference
0240     // 2- check on lastTSOS-lastButOneUpdatedTSOS eta difference
0241     // After this choice:
0242     // A- extract the measurements compatible with the initialTSOS (seed)
0243     // B- extract the measurements compatible with the lastButOneUpdatedTSOS
0244     // In ORCA the choice was 1A. Here I will try 1B and if it fail I'll try 1A
0245     // another possibility could be 2B and then 1A.
0246 
0247     // if no measurement found and the current TSOS has an eta very different
0248     // wrt the initial one (i.e. seed), then try to find the measurements
0249     // according to the lastButOne FTS. (1B)
0250     double lastdEta = fabs(lastTSOS.freeTrajectoryState()->momentum().eta() - eta0);
0251     if (bestMeasurements.empty() && lastdEta > 0.1) {
0252       LogTrace(metname) << "No measurement and big eta variation wrt seed" << endl
0253                         << "trying with lastButOneUpdatedTSOS";
0254       bestMeasurements = findBestMeasurements(*layer, theLastButOneUpdatedTSOS);
0255     }
0256 
0257     //if no measurement found and the current FTS has an eta very different
0258     //wrt the initial one (i.e. seed), then try to find the measurements
0259     //according to the initial FTS. (1A)
0260     if (bestMeasurements.empty() && lastdEta > 0.1) {
0261       LogTrace(metname) << "No measurement and big eta variation wrt seed" << endl << "tryng with seed TSOS";
0262       bestMeasurements = findBestMeasurements(*layer, initialTSOS);
0263     }
0264 
0265     // FIXME: uncomment this line!!
0266     // if(!bestMeasurement && firstTime) break;
0267 
0268     if (!bestMeasurements.empty()) {
0269       incrementCompatibleChamberCounters(*layer);
0270       bool added = false;
0271       for (std::vector<TrajectoryMeasurement>::const_iterator tmItr = bestMeasurements.begin();
0272            tmItr != bestMeasurements.end();
0273            ++tmItr) {
0274         added |= update(*layer, &(*tmItr), trajectory);
0275         lastTSOS = theLastUpdatedTSOS;
0276       }
0277       if (added) {
0278         incrementChamberCounters(*layer);
0279         theDetLayers.push_back(*layer);
0280       }
0281     }
0282     // SL in case no valid mesurement is found, still I want to use the predicted
0283     // state for the following measurement serches. I take the first in the
0284     // container. FIXME!!! I want to carefully check this!!!!!
0285     else {
0286       LogTrace(metname) << "No best measurement found" << endl;
0287       //       if (!theMeasurementCache.empty()){
0288       //    LogTrace(metname)<<"but the #of measurement is "<<theMeasurementCache.size()<<endl;
0289       //         lastTSOS = theMeasurementCache.front().predictedState();
0290       //       }
0291     }
0292   }  // loop over layers
0293 }
0294 
0295 std::vector<TrajectoryMeasurement> StandAloneMuonFilter::findBestMeasurements(const DetLayer* layer,
0296                                                                               const TrajectoryStateOnSurface& tsos) {
0297   const std::string metname = "Muon|RecoMuon|StandAloneMuonFilter";
0298 
0299   std::vector<TrajectoryMeasurement> result;
0300   std::vector<TrajectoryMeasurement> measurements;
0301 
0302   if (theOverlappingChambersFlag && layer->hasGroups()) {
0303     std::vector<TrajectoryMeasurementGroup> measurementGroups =
0304         theMeasurementExtractor->groupedMeasurements(layer, tsos, *propagator(), *estimator());
0305 
0306     if (theFitDirection == outsideIn) {
0307       LogTrace(metname) << "Reversing the order of groupedMeasurements as the direction of the fit is outside-in";
0308       reverse(measurementGroups.begin(), measurementGroups.end());
0309       // this should be fixed either in RecoMuon/MeasurementDet/MuonDetLayerMeasurements or
0310       // RecoMuon/DetLayers/MuRingForwardDoubleLayer
0311     }
0312 
0313     for (std::vector<TrajectoryMeasurementGroup>::const_iterator tmGroupItr = measurementGroups.begin();
0314          tmGroupItr != measurementGroups.end();
0315          ++tmGroupItr) {
0316       measurements = tmGroupItr->measurements();
0317       LogTrace(metname) << "Number of Trajectory Measurement: " << measurements.size();
0318 
0319       const TrajectoryMeasurement* bestMeasurement =
0320           bestMeasurementFinder()->findBestMeasurement(measurements, propagator());
0321 
0322       if (bestMeasurement)
0323         result.push_back(*bestMeasurement);
0324     }
0325   } else {
0326     measurements = theMeasurementExtractor->measurements(layer, tsos, *propagator(), *estimator());
0327     LogTrace(metname) << "Number of Trajectory Measurement: " << measurements.size();
0328     const TrajectoryMeasurement* bestMeasurement =
0329         bestMeasurementFinder()->findBestMeasurement(measurements, propagator());
0330     if (bestMeasurement)
0331       result.push_back(*bestMeasurement);
0332   }
0333   return result;
0334 }
0335 
0336 bool StandAloneMuonFilter::update(const DetLayer* layer, const TrajectoryMeasurement* meas, Trajectory& trajectory) {
0337   const std::string metname = "Muon|RecoMuon|StandAloneMuonFilter";
0338   MuonPatternRecoDumper debug;
0339 
0340   LogTrace(metname) << "best measurement found"
0341                     << "\n"
0342                     << "updating the trajectory..." << endl;
0343   pair<bool, TrajectoryStateOnSurface> result = updator()->update(meas, trajectory, propagator());
0344   LogTrace(metname) << "trajectory updated: " << result.first << endl;
0345   LogTrace(metname) << debug.dumpTSOS(result.second);
0346 
0347   if (result.first) {
0348     theLastButOneUpdatedTSOS = theLastUpdatedTSOS;
0349     theLastUpdatedTSOS = result.second;
0350   }
0351 
0352   if (result.second.isValid())
0353     theLastCompatibleTSOS = result.second;
0354 
0355   return result.first;
0356 }
0357 
0358 void StandAloneMuonFilter::createDefaultTrajectory(const Trajectory& oldTraj, Trajectory& defTraj) {
0359   Trajectory::DataContainer const& oldMeas = oldTraj.measurements();
0360   defTraj.reserve(oldMeas.size());
0361 
0362   for (Trajectory::DataContainer::const_iterator itm = oldMeas.begin(); itm != oldMeas.end(); itm++) {
0363     if (!(*itm).recHit()->isValid())
0364       defTraj.push(*itm, (*itm).estimate());
0365     else {
0366       MuonTransientTrackingRecHit::MuonRecHitPointer invRhPtr =
0367           MuonTransientTrackingRecHit::specificBuild((*itm).recHit()->det(), (*itm).recHit()->hit());
0368       invRhPtr->invalidateHit();
0369       TrajectoryMeasurement invRhMeas(
0370           (*itm).forwardPredictedState(), (*itm).updatedState(), invRhPtr, (*itm).estimate(), (*itm).layer());
0371       defTraj.push(std::move(invRhMeas), (*itm).estimate());
0372     }
0373 
0374   }  // end for
0375 }