MeasurementDet/src/MTDDetLayerMeasurements.cc

0001 /** \class MTDDetLayerMeasurements
0002  *  The class to access recHits and TrajectoryMeasurements from DetLayer.
0003  *
0004  *  \author B. Tannenwald
0005  *  Adapted from RecoMuon version.
0006  *
0007  */
0008
0009 #include "RecoMTD/MeasurementDet/interface/MTDDetLayerMeasurements.h"
0010
0011 #include "TrackingTools/PatternTools/interface/TrajectoryMeasurement.h"
0012 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0013
0014 #include "TrackingTools/PatternTools/interface/TrajMeasLessEstim.h"
0015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0016 #include "FWCore/ServiceRegistry/interface/Service.h"
0017
0018 typedef std::shared_ptr<GenericTransientTrackingRecHit> MTDRecHitPointer;
0019 typedef std::vector<GenericTransientTrackingRecHit::RecHitPointer> MTDRecHitContainer;
0020 typedef MTDDetLayerMeasurements::MeasurementContainer MeasurementContainer;
0021
0022 MTDDetLayerMeasurements::MTDDetLayerMeasurements(const edm::InputTag& mtdlabel, edm::ConsumesCollector& iC)
0023     : theMTDToken(iC.consumes<MTDTrackingRecHit>(mtdlabel)),
0024       theMTDRecHits(),
0025       theMTDEventCacheID(0),
0026       theEvent(nullptr) {}
0027
0028 MTDDetLayerMeasurements::~MTDDetLayerMeasurements() {}
0029
0030 MTDRecHitContainer MTDDetLayerMeasurements::recHits(const GeomDet* geomDet, const edm::Event& iEvent) {
0031   DetId geoId = geomDet->geographicalId();
0032   theEvent = &iEvent;
0033   MTDRecHitContainer result;
0034
0035   checkMTDRecHits();
0036
0037   // Create the ChamberId
0038   DetId detId(geoId.rawId());
0039   LogDebug("MTDDetLayerMeasurements") << "(MTD): " << static_cast<MTDDetId>(detId) << std::endl;
0040
0041   // Get the MTD-Segment which relies on this chamber
0042   auto detset = (*theMTDRecHits)[detId];
0043
0044   for (const auto& rechit : detset)
0045     result.push_back(GenericTransientTrackingRecHit::build(geomDet, &rechit));
0046
0047   return result;
0048 }
0049
0050 void MTDDetLayerMeasurements::checkMTDRecHits() {
0051   LogDebug("MTDDetLayerMeasurements") << "Checking MTD RecHits";
0052   checkEvent();
0053   auto cacheID = theEvent->cacheIdentifier();
0054   if (cacheID == theMTDEventCacheID)
0055     return;
0056
0057   {
0058     theEvent->getByToken(theMTDToken, theMTDRecHits);
0059     theMTDEventCacheID = cacheID;
0060   }
0061   if (!theMTDRecHits.isValid()) {
0062     throw cms::Exception("MTDDetLayerMeasurements") << "Cannot get MTD RecHits";
0063   }
0064 }
0065
0066 template <class T>
0067 T MTDDetLayerMeasurements::sortResult(T& result) {
0068   if (!result.empty()) {
0069     sort(result.begin(), result.end(), TrajMeasLessEstim());
0070   }
0071
0072   return result;
0073 }
0074
0075 ///measurements method if already got the Event
0076 MeasurementContainer MTDDetLayerMeasurements::measurements(const DetLayer* layer,
0077                                                            const TrajectoryStateOnSurface& startingState,
0078                                                            const Propagator& prop,
0079                                                            const MeasurementEstimator& est) {
0080   checkEvent();
0081   return measurements(layer, startingState, prop, est, *theEvent);
0082 }
0083
0084 MeasurementContainer MTDDetLayerMeasurements::measurements(const DetLayer* layer,
0085                                                            const TrajectoryStateOnSurface& startingState,
0086                                                            const Propagator& prop,
0087                                                            const MeasurementEstimator& est,
0088                                                            const edm::Event& iEvent) {
0089   MeasurementContainer result;
0090
0091   const auto& dss = layer->compatibleDets(startingState, prop, est);
0092   LogDebug("MTDDetLayerMeasurements") << "compatibleDets: " << dss.size() << std::endl;
0093
0094   for (const auto& dws : dss) {
0095     MeasurementContainer detMeasurements = measurements(layer, dws.first, dws.second, est, iEvent);
0096     result.insert(result.end(), detMeasurements.begin(), detMeasurements.end());
0097   }
0098
0099   return sortResult(result);
0100 }
0101
0102 MeasurementContainer MTDDetLayerMeasurements::measurements(const DetLayer* layer,
0103                                                            const GeomDet* det,
0104                                                            const TrajectoryStateOnSurface& stateOnDet,
0105                                                            const MeasurementEstimator& est,
0106                                                            const edm::Event& iEvent) {
0107   MeasurementContainer result;
0108
0109   // Get the Segments which relies on the GeomDet given by compatibleDets
0110   MTDRecHitContainer mtdRecHits = recHits(det, iEvent);
0111
0112   // Create the Trajectory Measurement
0113   for (const auto& rechit : mtdRecHits) {
0114     MeasurementEstimator::HitReturnType estimate = est.estimate(stateOnDet, *rechit);
0115     LogDebug("RecoMTD") << "Dimension: " << rechit->dimension() << " Chi2: " << estimate.second << std::endl;
0116     if (estimate.first) {
0117       result.push_back(TrajectoryMeasurement(stateOnDet, rechit, estimate.second, layer));
0118     }
0119   }
0120
0121   return sortResult(result);
0122 }
0123
0124 MeasurementContainer MTDDetLayerMeasurements::fastMeasurements(const DetLayer* layer,
0125                                                                const TrajectoryStateOnSurface& theStateOnDet,
0126                                                                const TrajectoryStateOnSurface& startingState,
0127                                                                const Propagator& prop,
0128                                                                const MeasurementEstimator& est,
0129                                                                const edm::Event& iEvent) {
0130   MeasurementContainer result;
0131   MTDRecHitContainer rhs = recHits(layer, iEvent);
0132   for (const auto& irh : rhs) {
0133     MeasurementEstimator::HitReturnType estimate = est.estimate(theStateOnDet, (*irh));
0134     if (estimate.first) {
0135       result.push_back(TrajectoryMeasurement(theStateOnDet, irh, estimate.second, layer));
0136     }
0137   }
0138
0139   return sortResult(result);
0140 }
0141
0142 ///fastMeasurements method if already got the Event
0143 MeasurementContainer MTDDetLayerMeasurements::fastMeasurements(const DetLayer* layer,
0144                                                                const TrajectoryStateOnSurface& theStateOnDet,
0145                                                                const TrajectoryStateOnSurface& startingState,
0146                                                                const Propagator& prop,
0147                                                                const MeasurementEstimator& est) {
0148   checkEvent();
0149   return fastMeasurements(layer, theStateOnDet, startingState, prop, est, *theEvent);
0150 }
0151
0152 std::vector<TrajectoryMeasurementGroup> MTDDetLayerMeasurements::groupedMeasurements(
0153     const DetLayer* layer,
0154     const TrajectoryStateOnSurface& startingState,
0155     const Propagator& prop,
0156     const MeasurementEstimator& est) {
0157   checkEvent();
0158   return groupedMeasurements(layer, startingState, prop, est, *theEvent);
0159 }
0160
0161 std::vector<TrajectoryMeasurementGroup> MTDDetLayerMeasurements::groupedMeasurements(
0162     const DetLayer* layer,
0163     const TrajectoryStateOnSurface& startingState,
0164     const Propagator& prop,
0165     const MeasurementEstimator& est,
0166     const edm::Event& iEvent) {
0167   std::vector<TrajectoryMeasurementGroup> result;
0168   // if we want to use the concept of InvalidRecHits,
0169   // we can reuse LayerMeasurements from TrackingTools/MeasurementDet
0170   std::vector<DetGroup> groups(layer->groupedCompatibleDets(startingState, prop, est));
0171
0172   for (const auto& grp : groups) {
0173     std::vector<TrajectoryMeasurement> groupMeasurements;
0174     for (const auto& detAndStateItr : grp) {
0175       std::vector<TrajectoryMeasurement> detMeasurements =
0176           measurements(layer, detAndStateItr.det(), detAndStateItr.trajectoryState(), est, iEvent);
0177       groupMeasurements.insert(groupMeasurements.end(), detMeasurements.begin(), detMeasurements.end());
0178     }
0179
0180     result.push_back(TrajectoryMeasurementGroup(sortResult(groupMeasurements), grp));
0181   }
0182
0183   return result;
0184 }
0185
0186 ///set event
0187 void MTDDetLayerMeasurements::setEvent(const edm::Event& event) { theEvent = &event; }
0188
0189 void MTDDetLayerMeasurements::checkEvent() const {
0190   if (!theEvent)
0191     throw cms::Exception("MTDDetLayerMeasurements") << "The event has not been set";
0192 }
0193
0194 MTDRecHitContainer MTDDetLayerMeasurements::recHits(const DetLayer* layer, const edm::Event& iEvent) {
0195   MTDRecHitContainer rhs;
0196
0197   std::vector<const GeomDet*> gds = layer->basicComponents();
0198
0199   for (const GeomDet* igd : gds) {
0200     MTDRecHitContainer detHits = recHits(igd, iEvent);
0201     rhs.insert(rhs.end(), detHits.begin(), detHits.end());
0202   }
0203   return rhs;
0204 }
0205
0206 MTDRecHitContainer MTDDetLayerMeasurements::recHits(const DetLayer* layer) {
0207   checkEvent();
0208   return recHits(layer, *theEvent);
0209 }