MeasurementDet/src/MuonDetLayerMeasurements.cc

0001 /** \class MuonDetLayerMeasurements
0002  *  The class to access recHits and TrajectoryMeasurements from DetLayer.
0003  *
0004  *  \author C. Liu, R. Bellan, N. Amapane
0005  *   \modified by C. Calabria to include GEMs
0006  *  \modified by D. Nash to include ME0s
0007  *
0008  *  \modified by C. Calabria & A.Sharma to include GEMs
0009  *
0010  *
0011  */
0012
0013 #include "RecoMuon/MeasurementDet/interface/MuonDetLayerMeasurements.h"
0014
0015 #include "TrackingTools/PatternTools/interface/TrajectoryMeasurement.h"
0016 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0017
0018 #include "TrackingTools/PatternTools/interface/TrajMeasLessEstim.h"
0019 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0020 #include "FWCore/ServiceRegistry/interface/Service.h"
0021
0022 typedef MuonTransientTrackingRecHit::MuonRecHitPointer MuonRecHitPointer;
0023 typedef MuonTransientTrackingRecHit::MuonRecHitContainer MuonRecHitContainer;
0024
0025 MuonDetLayerMeasurements::MuonDetLayerMeasurements(edm::InputTag dtlabel,
0026                                                    edm::InputTag csclabel,
0027                                                    edm::InputTag rpclabel,
0028                                                    edm::InputTag gemlabel,
0029                                                    edm::InputTag me0label,
0030                                                    edm::ConsumesCollector& iC,
0031                                                    bool enableDT,
0032                                                    bool enableCSC,
0033                                                    bool enableRPC,
0034                                                    bool enableGEM,
0035                                                    bool enableME0)
0036     : enableDTMeasurement(enableDT),
0037       enableCSCMeasurement(enableCSC),
0038       enableRPCMeasurement(enableRPC),
0039       enableGEMMeasurement(enableGEM),
0040       enableME0Measurement(enableME0),
0041       theDTRecHits(),
0042       theCSCRecHits(),
0043       theRPCRecHits(),
0044       theGEMRecHits(),
0045       theDTEventCacheID(0),
0046       theCSCEventCacheID(0),
0047       theRPCEventCacheID(0),
0048       theGEMEventCacheID(0),
0049       theME0EventCacheID(0),
0050       theEvent(nullptr) {
0051   dtToken_ = iC.consumes<DTRecSegment4DCollection>(dtlabel);
0052   cscToken_ = iC.consumes<CSCSegmentCollection>(csclabel);
0053   rpcToken_ = iC.consumes<RPCRecHitCollection>(rpclabel);
0054   gemToken_ = iC.consumes<GEMRecHitCollection>(gemlabel);
0055   me0Token_ = iC.consumes<ME0SegmentCollection>(me0label);
0056
0057   static std::atomic<int> procInstance{0};
0058   std::ostringstream sDT;
0059   sDT << "MuonDetLayerMeasurements::checkDTRecHits::" << procInstance;
0060   //  theDTCheckName = sDT.str();
0061   std::ostringstream sRPC;
0062   sRPC << "MuonDetLayerMeasurements::checkRPCRecHits::" << procInstance;
0063   //theRPCCheckName = sRPC.str();
0064   std::ostringstream sCSC;
0065   sCSC << "MuonDetLayerMeasurements::checkCSCRecHits::" << procInstance;
0066   //theCSCCheckName = sCSC.str();
0067   std::ostringstream sGEM;
0068   sGEM << "MuonDetLayerMeasurements::checkGEMRecHits::" << procInstance;
0069   //theGEMCheckName = sGEM.str();
0070   std::ostringstream sME0;
0071   sME0 << "MuonDetLayerMeasurements::checkME0RecHits::" << procInstance;
0072   //theME0CheckName = sME0.str();
0073   procInstance++;
0074 }
0075
0076 MuonDetLayerMeasurements::~MuonDetLayerMeasurements() {}
0077
0078 MuonRecHitContainer MuonDetLayerMeasurements::recHits(const GeomDet* geomDet, const edm::Event& iEvent) {
0079   DetId geoId = geomDet->geographicalId();
0080   theEvent = &iEvent;
0081   MuonRecHitContainer result;
0082
0083   if (geoId.subdetId() == MuonSubdetId::DT) {
0084     if (enableDTMeasurement) {
0085       checkDTRecHits();
0086
0087       // Create the ChamberId
0088       DTChamberId chamberId(geoId.rawId());
0089       LogDebug("Muon|RecoMuon|MuonDetLayerMeasurements") << "(DT): " << chamberId << std::endl;
0090
0091       // Get the DT-Segment which relies on this chamber
0092       DTRecSegment4DCollection::range range = theDTRecHits->get(chamberId);
0093
0094       // Create the MuonTransientTrackingRechit
0095       for (DTRecSegment4DCollection::const_iterator rechit = range.first; rechit != range.second; ++rechit)
0096         result.push_back(MuonTransientTrackingRecHit::specificBuild(geomDet, &*rechit));
0097     }
0098   }
0099
0100   else if (geoId.subdetId() == MuonSubdetId::CSC) {
0101     if (enableCSCMeasurement) {
0102       checkCSCRecHits();
0103
0104       // Create the chamber Id
0105       CSCDetId chamberId(geoId.rawId());
0106       LogDebug("Muon|RecoMuon|MuonDetLayerMeasurements") << "(CSC): " << chamberId << std::endl;
0107
0108       // Get the CSC-Segment which relies on this chamber
0109       CSCSegmentCollection::range range = theCSCRecHits->get(chamberId);
0110
0111       // Create the MuonTransientTrackingRecHit
0112       for (CSCSegmentCollection::const_iterator rechit = range.first; rechit != range.second; ++rechit)
0113         result.push_back(MuonTransientTrackingRecHit::specificBuild(geomDet, &*rechit));
0114     }
0115   }
0116
0117   else if (geoId.subdetId() == MuonSubdetId::RPC) {
0118     if (enableRPCMeasurement) {
0119       checkRPCRecHits();
0120
0121       // Create the chamber Id
0122       RPCDetId chamberId(geoId.rawId());
0123       LogDebug("Muon|RecoMuon|MuonDetLayerMeasurements") << "(RPC): " << chamberId << std::endl;
0124
0125       // Get the RPC-Segment which relies on this chamber
0126       RPCRecHitCollection::range range = theRPCRecHits->get(chamberId);
0127
0128       // Create the MuonTransientTrackingRecHit
0129       for (RPCRecHitCollection::const_iterator rechit = range.first; rechit != range.second; ++rechit)
0130         result.push_back(MuonTransientTrackingRecHit::specificBuild(geomDet, &*rechit));
0131     }
0132   } else if (geoId.subdetId() == MuonSubdetId::GEM) {
0133     if (enableGEMMeasurement) {
0134       checkGEMRecHits();
0135
0136       // Create the chamber Id
0137       GEMDetId chamberId(geoId.rawId());
0138
0139       LogDebug("Muon|RecoMuon|MuonDetLayerMeasurements") << "(GEM): " << chamberId << std::endl;
0140
0141       // Get the GEM-Segment which relies on this chamber
0142       GEMRecHitCollection::range range = theGEMRecHits->get(chamberId);
0143
0144       LogDebug("Muon|RecoMuon|MuonDetLayerMeasurements")
0145           << "Number of GEM rechits available =  " << theGEMRecHits->size() << ", from chamber: " << chamberId
0146           << std::endl;
0147
0148       // Create the MuonTransientTrackingRecHit
0149       for (GEMRecHitCollection::const_iterator rechit = range.first; rechit != range.second; ++rechit)
0150         result.push_back(MuonTransientTrackingRecHit::specificBuild(geomDet, &*rechit));
0151       LogDebug("Muon|RecoMuon|MuonDetLayerMeasurements") << "Number of GEM rechits = " << result.size() << std::endl;
0152     }
0153   }
0154
0155   else if (geoId.subdetId() == MuonSubdetId::ME0) {
0156     LogDebug("Muon|RecoMuon|MuonDetLayerMeasurements") << "(ME0): identified" << std::endl;
0157     if (enableME0Measurement) {
0158       checkME0RecHits();
0159
0160       // Create the chamber Id
0161       ME0DetId chamberId(geoId.rawId());
0162
0163       // Get the ME0-Segment which relies on this chamber
0164       // Getting rechits right now, not segments - maybe it should be segments?
0165       ME0SegmentCollection::range range = theME0RecHits->get(chamberId);
0166
0167       LogDebug("Muon|RecoMuon|MuonDetLayerMeasurements")
0168           << "Number of ME0 rechits available =  " << theME0RecHits->size() << ", from chamber: " << chamberId
0169           << std::endl;
0170
0171       // Create the MuonTransientTrackingRecHit
0172       for (ME0SegmentCollection::const_iterator rechit = range.first; rechit != range.second; ++rechit) {
0173         LogDebug("Muon|RecoMuon|MuonDetLayerMeasurements") << "On ME0 iteration " << std::endl;
0174         result.push_back(MuonTransientTrackingRecHit::specificBuild(geomDet, &*rechit));
0175       }
0176       LogDebug("Muon|RecoMuon|MuonDetLayerMeasurements") << "Number of ME0 rechits = " << result.size() << std::endl;
0177     }
0178   } else {
0179     // wrong type
0180     throw cms::Exception("MuonDetLayerMeasurements") << "The DetLayer with det " << geoId.det() << " subdet "
0181                                                      << geoId.subdetId() << " is not a valid Muon DetLayer. ";
0182     LogDebug("Muon|RecoMuon|MuonDetLayerMeasurements") << "The DetLayer with det " << geoId.det() << " subdet "
0183                                                        << geoId.subdetId() << " is not a valid Muon DetLayer. ";
0184   }
0185   if (enableME0Measurement) {
0186     LogDebug("Muon|RecoMuon|MuonDetLayerMeasurements") << "(ME0): enabled" << std::endl;
0187   }
0188
0189   if (enableGEMMeasurement) {
0190     LogDebug("Muon|RecoMuon|MuonDetLayerMeasurements") << "(GEM): enabled" << std::endl;
0191   }
0192   return result;
0193 }
0194
0195 void MuonDetLayerMeasurements::checkDTRecHits() {
0196   checkEvent();
0197   auto const cacheID = theEvent->cacheIdentifier();
0198   if (cacheID == theDTEventCacheID)
0199     return;
0200
0201   { theEvent->getByToken(dtToken_, theDTRecHits); }
0202   if (!theDTRecHits.isValid()) {
0203     throw cms::Exception("MuonDetLayerMeasurements") << "Cannot get DT RecHits";
0204   }
0205 }
0206
0207 void MuonDetLayerMeasurements::checkCSCRecHits() {
0208   checkEvent();
0209   auto cacheID = theEvent->cacheIdentifier();
0210   if (cacheID == theCSCEventCacheID)
0211     return;
0212
0213   {
0214     theEvent->getByToken(cscToken_, theCSCRecHits);
0215     theCSCEventCacheID = cacheID;
0216   }
0217   if (!theCSCRecHits.isValid()) {
0218     throw cms::Exception("MuonDetLayerMeasurements") << "Cannot get CSC RecHits";
0219   }
0220 }
0221
0222 void MuonDetLayerMeasurements::checkRPCRecHits() {
0223   checkEvent();
0224   auto cacheID = theEvent->cacheIdentifier();
0225   if (cacheID == theRPCEventCacheID)
0226     return;
0227
0228   {
0229     theEvent->getByToken(rpcToken_, theRPCRecHits);
0230     theRPCEventCacheID = cacheID;
0231   }
0232   if (!theRPCRecHits.isValid()) {
0233     throw cms::Exception("MuonDetLayerMeasurements") << "Cannot get RPC RecHits";
0234   }
0235 }
0236
0237 void MuonDetLayerMeasurements::checkGEMRecHits() {
0238   checkEvent();
0239   auto cacheID = theEvent->cacheIdentifier();
0240   if (cacheID == theGEMEventCacheID)
0241     return;
0242
0243   {
0244     theEvent->getByToken(gemToken_, theGEMRecHits);
0245     theGEMEventCacheID = cacheID;
0246   }
0247   if (!theGEMRecHits.isValid()) {
0248     throw cms::Exception("MuonDetLayerMeasurements") << "Cannot get GEM RecHits";
0249   }
0250 }
0251
0252 void MuonDetLayerMeasurements::checkME0RecHits() {
0253   LogDebug("Muon|RecoMuon|MuonDetLayerMeasurements") << "Checking ME0 RecHits";
0254   checkEvent();
0255   auto cacheID = theEvent->cacheIdentifier();
0256   if (cacheID == theME0EventCacheID)
0257     return;
0258
0259   {
0260     theEvent->getByToken(me0Token_, theME0RecHits);
0261     theME0EventCacheID = cacheID;
0262   }
0263   if (!theME0RecHits.isValid()) {
0264     throw cms::Exception("MuonDetLayerMeasurements") << "Cannot get ME0 RecHits";
0265     LogDebug("Muon|RecoMuon|MuonDetLayerMeasurements") << "Cannot get ME0 RecHits";
0266   }
0267 }
0268
0269 ///measurements method if already got the Event
0270 MeasurementContainer MuonDetLayerMeasurements::measurements(const DetLayer* layer,
0271                                                             const TrajectoryStateOnSurface& startingState,
0272                                                             const Propagator& prop,
0273                                                             const MeasurementEstimator& est) {
0274   checkEvent();
0275   return measurements(layer, startingState, prop, est, *theEvent);
0276 }
0277
0278 MeasurementContainer MuonDetLayerMeasurements::measurements(const DetLayer* layer,
0279                                                             const TrajectoryStateOnSurface& startingState,
0280                                                             const Propagator& prop,
0281                                                             const MeasurementEstimator& est,
0282                                                             const edm::Event& iEvent) {
0283   MeasurementContainer result;
0284
0285   std::vector<DetWithState> dss = layer->compatibleDets(startingState, prop, est);
0286   LogDebug("RecoMuon") << "compatibleDets: " << dss.size() << std::endl;
0287
0288   for (std::vector<DetWithState>::const_iterator detWithStateItr = dss.begin(); detWithStateItr != dss.end();
0289        ++detWithStateItr) {
0290     MeasurementContainer detMeasurements =
0291         measurements(layer, detWithStateItr->first, detWithStateItr->second, est, iEvent);
0292     result.insert(result.end(), detMeasurements.begin(), detMeasurements.end());
0293   }
0294
0295   if (!result.empty())
0296     sort(result.begin(), result.end(), TrajMeasLessEstim());
0297
0298   return result;
0299 }
0300
0301 MeasurementContainer MuonDetLayerMeasurements::measurements(const DetLayer* layer,
0302                                                             const GeomDet* det,
0303                                                             const TrajectoryStateOnSurface& stateOnDet,
0304                                                             const MeasurementEstimator& est,
0305                                                             const edm::Event& iEvent) {
0306   MeasurementContainer result;
0307
0308   DetId geoId = det->geographicalId();
0309   //no chamber here that is actually an me0....
0310   if (geoId.subdetId() == MuonSubdetId::ME0) {
0311     if (enableME0Measurement) {
0312       ME0DetId chamberId(geoId.rawId());
0313       LogDebug("Muon|RecoMuon|MuonDetLayerMeasurements")
0314           << "ME0 Chamber ID in measurements: " << chamberId << std::endl;
0315     }
0316   }
0317
0318   // Get the Segments which relies on the GeomDet given by compatibleDets
0319   MuonRecHitContainer muonRecHits = recHits(det, iEvent);
0320
0321   // Create the Trajectory Measurement
0322   for (MuonRecHitContainer::const_iterator rechit = muonRecHits.begin(); rechit != muonRecHits.end(); ++rechit) {
0323     MeasurementEstimator::HitReturnType estimate = est.estimate(stateOnDet, **rechit);
0324     LogDebug("RecoMuon") << "Dimension: " << (*rechit)->dimension() << " Chi2: " << estimate.second << std::endl;
0325     if (estimate.first) {
0326       result.push_back(TrajectoryMeasurement(stateOnDet, *rechit, estimate.second, layer));
0327     }
0328   }
0329
0330   if (!result.empty())
0331     sort(result.begin(), result.end(), TrajMeasLessEstim());
0332
0333   return result;
0334 }
0335
0336 MeasurementContainer MuonDetLayerMeasurements::fastMeasurements(const DetLayer* layer,
0337                                                                 const TrajectoryStateOnSurface& theStateOnDet,
0338                                                                 const TrajectoryStateOnSurface& startingState,
0339                                                                 const Propagator& prop,
0340                                                                 const MeasurementEstimator& est,
0341                                                                 const edm::Event& iEvent) {
0342   MeasurementContainer result;
0343   MuonRecHitContainer rhs = recHits(layer, iEvent);
0344   for (MuonRecHitContainer::const_iterator irh = rhs.begin(); irh != rhs.end(); irh++) {
0345     MeasurementEstimator::HitReturnType estimate = est.estimate(theStateOnDet, (**irh));
0346     if (estimate.first) {
0347       result.push_back(TrajectoryMeasurement(theStateOnDet, (*irh), estimate.second, layer));
0348     }
0349   }
0350
0351   if (!result.empty()) {
0352     sort(result.begin(), result.end(), TrajMeasLessEstim());
0353   }
0354
0355   return result;
0356 }
0357
0358 ///fastMeasurements method if already got the Event
0359 MeasurementContainer MuonDetLayerMeasurements::fastMeasurements(const DetLayer* layer,
0360                                                                 const TrajectoryStateOnSurface& theStateOnDet,
0361                                                                 const TrajectoryStateOnSurface& startingState,
0362                                                                 const Propagator& prop,
0363                                                                 const MeasurementEstimator& est) {
0364   checkEvent();
0365   return fastMeasurements(layer, theStateOnDet, startingState, prop, est, *theEvent);
0366 }
0367
0368 std::vector<TrajectoryMeasurementGroup> MuonDetLayerMeasurements::groupedMeasurements(
0369     const DetLayer* layer,
0370     const TrajectoryStateOnSurface& startingState,
0371     const Propagator& prop,
0372     const MeasurementEstimator& est) {
0373   checkEvent();
0374   return groupedMeasurements(layer, startingState, prop, est, *theEvent);
0375 }
0376
0377 std::vector<TrajectoryMeasurementGroup> MuonDetLayerMeasurements::groupedMeasurements(
0378     const DetLayer* layer,
0379     const TrajectoryStateOnSurface& startingState,
0380     const Propagator& prop,
0381     const MeasurementEstimator& est,
0382     const edm::Event& iEvent) {
0383   std::vector<TrajectoryMeasurementGroup> result;
0384   // if we want to use the concept of InvalidRecHits,
0385   // we can reuse LayerMeasurements from TrackingTools/MeasurementDet
0386   std::vector<DetGroup> groups(layer->groupedCompatibleDets(startingState, prop, est));
0387
0388   // this should be fixed either in RecoMuon/MeasurementDet/MuonDetLayerMeasurements or
0389   // RecoMuon/DetLayers/MuRingForwardDoubleLayer
0390   // and removed the reverse operation in StandAloneMuonFilter::findBestMeasurements
0391
0392   for (std::vector<DetGroup>::const_iterator grp = groups.begin(); grp != groups.end(); ++grp) {
0393     std::vector<TrajectoryMeasurement> groupMeasurements;
0394     for (DetGroup::const_iterator detAndStateItr = grp->begin(); detAndStateItr != grp->end(); ++detAndStateItr) {
0395       std::vector<TrajectoryMeasurement> detMeasurements =
0396           measurements(layer, detAndStateItr->det(), detAndStateItr->trajectoryState(), est, iEvent);
0397       groupMeasurements.insert(groupMeasurements.end(), detMeasurements.begin(), detMeasurements.end());
0398     }
0399
0400     if (!groupMeasurements.empty())
0401       std::sort(groupMeasurements.begin(), groupMeasurements.end(), TrajMeasLessEstim());
0402
0403     result.push_back(TrajectoryMeasurementGroup(groupMeasurements, *grp));
0404   }
0405
0406   return result;
0407 }
0408
0409 ///set event
0410 void MuonDetLayerMeasurements::setEvent(const edm::Event& event) { theEvent = &event; }
0411
0412 void MuonDetLayerMeasurements::checkEvent() const {
0413   if (!theEvent)
0414     throw cms::Exception("MuonDetLayerMeasurements") << "The event has not been set";
0415 }
0416
0417 MuonRecHitContainer MuonDetLayerMeasurements::recHits(const DetLayer* layer, const edm::Event& iEvent) {
0418   MuonRecHitContainer rhs;
0419
0420   std::vector<const GeomDet*> gds = layer->basicComponents();
0421
0422   for (std::vector<const GeomDet*>::const_iterator igd = gds.begin(); igd != gds.end(); igd++) {
0423     MuonRecHitContainer detHits = recHits(*igd, iEvent);
0424     rhs.insert(rhs.end(), detHits.begin(), detHits.end());
0425   }
0426   return rhs;
0427 }
0428
0429 MuonRecHitContainer MuonDetLayerMeasurements::recHits(const DetLayer* layer) {
0430   checkEvent();
0431   return recHits(layer, *theEvent);
0432 }