MeasurementDet/src/LayerMeasurements.cc

0001 #include "TrackingTools/MeasurementDet/interface/LayerMeasurements.h"
0002 #include "TrackingTools/PatternTools/interface/TrajectoryMeasurement.h"
0003 #include "TrackingTools/PatternTools/interface/TrajMeasLessEstim.h"
0004 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0005
0006 #include "TrackingTools/MeasurementDet/interface/MeasurementDetException.h"
0007 #include "TrackingTools/MeasurementDet/interface/MeasurementDetSystem.h"
0008 #include "TrackingTools/MeasurementDet/interface/MeasurementDet.h"
0009 #include "TrackingTools/MeasurementDet/interface/TrajectoryMeasurementGroup.h"
0010
0011 #include "TrackingTools/DetLayers/interface/GeometricSearchDet.h"
0012 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0013 #include "TrackingTools/DetLayers/interface/DetGroup.h"
0014
0015 #include "TrackingTools/TransientTrackingRecHit/interface/InvalidTransientRecHit.h"
0016
0017 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0018 #include "FWCore/Utilities/interface/Likely.h"
0019
0020 #include <algorithm>
0021 #include <cstdlib>
0022
0023 using namespace std;
0024
0025 namespace {
0026   typedef GeometricSearchDet::DetWithState DetWithState;
0027   inline void addInvalidMeas(std::vector<TrajectoryMeasurement>& result,
0028                              const TrajectoryStateOnSurface& ts,
0029                              const GeomDet& det,
0030                              const DetLayer& layer) {
0031     result.emplace_back(ts, std::make_shared<InvalidTrackingRecHit>(det, TrackingRecHit::missing), 0.F, &layer);
0032   }
0033
0034   // layerIsCompatibleWithTSOS:
0035   //  - enforce minimum requirements of compatibility between DetLayer and TrajectoryStateOnSurface
0036   //    (e.g. global position of TSOS within detector volume)
0037   //  - current implementation does not make use of layer argument,
0038   //    see below for an example of how to use it
0039   //      #include "Geometry/CommonDetUnit/interface/GeomDetEnumerators.h"
0040   //       [..]
0041   //      using namespace GeomDetEnumerators;
0042   //      auto const& det = layer.subDetector();
0043   //      if (det == PixelBarrel or det == PixelEndcap) {
0044   //        return xa < 2e1 and ya < 2e1 and za < 6e1;
0045   //      }
0046   bool layerIsCompatibleWithTSOS(DetLayer const& /* layer */, TrajectoryStateOnSurface const& tsos) {
0047     auto const& gp = tsos.globalPosition();
0048     auto const xa = std::abs(gp.x());
0049     auto const ya = std::abs(gp.y());
0050     auto const za = std::abs(gp.z());
0051     return xa < 1e3 and ya < 1e3 and za < 12e2;
0052   }
0053
0054   /** The std::vector<DetWithState> passed to this method should not be empty.
0055    *  In case of no compatible dets the result should be either an empty container if
0056    *  the det is itself incompatible, or a container with one invalid measurement
0057    *  on the det surface. The method does not have enough information to do
0058    *  this efficiently, so it should be done by the caller, or an exception will
0059    *  be thrown (DetLogicError).
0060    */
0061   inline std::vector<TrajectoryMeasurement> get(MeasurementDetSystem const& detSystem,
0062                                                 MeasurementTrackerEvent const& data,
0063                                                 const DetLayer& layer,
0064                                                 std::vector<DetWithState> const& compatDets,
0065                                                 const TrajectoryStateOnSurface& ts,
0066                                                 const Propagator& prop,
0067                                                 const MeasurementEstimator& est) {
0068     std::vector<TrajectoryMeasurement> result;
0069     typedef TrajectoryMeasurement TM;
0070
0071     tracking::TempMeasurements tmps;
0072
0073     for (auto const& ds : compatDets) {
0074       if (not layerIsCompatibleWithTSOS(layer, ds.second))
0075         continue;
0076
0077       MeasurementDetWithData mdet = detSystem.idToDet(ds.first->geographicalId(), data);
0078       if UNLIKELY (mdet.isNull()) {
0079         throw MeasurementDetException("MeasurementDet not found");
0080       }
0081
0082       if (mdet.measurements(ds.second, est, tmps))
0083         for (std::size_t i = 0; i != tmps.size(); ++i)
0084           result.emplace_back(ds.second, std::move(tmps.hits[i]), tmps.distances[i], &layer);
0085       tmps.clear();
0086     }
0087     // WARNING: we might end up with more than one invalid hit of type 'inactive' in result
0088     // to be fixed in order to avoid usless double traj candidates.
0089
0090     // sort the final result
0091     if (result.size() > 1) {
0092       sort(result.begin(), result.end(), TrajMeasLessEstim());
0093     }
0094
0095     if (!result.empty()) {
0096       // invalidMeas on Det of most compatible hit
0097       addInvalidMeas(result, result.front().predictedState(), *(result.front().recHit()->det()), layer);
0098     } else {
0099       // invalid state on first compatible Det
0100       addInvalidMeas(result, compatDets.front().second, *(compatDets.front().first), layer);
0101     }
0102
0103     return result;
0104   }
0105
0106   void addInvalidMeas(vector<TrajectoryMeasurement>& measVec, const DetGroup& group, const DetLayer& layer) {
0107     if (!measVec.empty()) {
0108       // invalidMeas on Det of most compatible hit
0109       auto const& ts = measVec.front().predictedState();
0110       auto toll = measVec.front().recHitR().det()->surface().bounds().significanceInside(
0111           ts.localPosition(), ts.localError().positionError());
0112       measVec.emplace_back(
0113           measVec.front().predictedState(),
0114           std::make_shared<InvalidTrackingRecHit>(*measVec.front().recHitR().det(), TrackingRecHit::missing),
0115           toll,
0116           &layer);
0117     } else if (!group.empty()) {
0118       // invalid state on first compatible Det
0119       auto const& ts = group.front().trajectoryState();
0120       auto toll = group.front().det()->surface().bounds().significanceInside(ts.localPosition(),
0121                                                                              ts.localError().positionError());
0122       measVec.emplace_back(group.front().trajectoryState(),
0123                            std::make_shared<InvalidTrackingRecHit>(*group.front().det(), TrackingRecHit::missing),
0124                            toll,
0125                            &layer);
0126     }
0127   }
0128
0129 }  // namespace
0130
0131 // return just valid hits, no sorting (for seeding mostly)
0132 std::vector<BaseTrackerRecHit*> LayerMeasurements::recHits(const DetLayer& layer,
0133                                                            const TrajectoryStateOnSurface& startingState,
0134                                                            const Propagator& prop,
0135                                                            const MeasurementEstimator& est) const {
0136   std::vector<BaseTrackerRecHit*> result;
0137   auto const& compatDets = layer.compatibleDets(startingState, prop, est);
0138   if (compatDets.empty())
0139     return result;
0140   for (auto const& ds : compatDets) {
0141     if (not layerIsCompatibleWithTSOS(layer, ds.second))
0142       continue;
0143
0144     auto mdet = detSystem_.idToDet(ds.first->geographicalId(), data_);
0145     mdet.recHits(result, ds.second, est);
0146   }
0147   return result;
0148 }
0149
0150 vector<TrajectoryMeasurement> LayerMeasurements::measurements(const DetLayer& layer,
0151                                                               const TrajectoryStateOnSurface& startingState,
0152                                                               const Propagator& prop,
0153                                                               const MeasurementEstimator& est) const {
0154   typedef DetLayer::DetWithState DetWithState;
0155
0156   vector<DetWithState> const& compatDets = layer.compatibleDets(startingState, prop, est);
0157
0158   if (!compatDets.empty())
0159     return get(detSystem_, data_, layer, compatDets, startingState, prop, est);
0160
0161   vector<TrajectoryMeasurement> result;
0162   pair<bool, TrajectoryStateOnSurface> compat = layer.compatible(startingState, prop, est);
0163
0164   if (compat.first and layerIsCompatibleWithTSOS(layer, compat.second)) {
0165     result.push_back(
0166         TrajectoryMeasurement(compat.second,
0167                               std::make_shared<InvalidTrackingRecHitNoDet>(layer.surface(), TrackingRecHit::inactive),
0168                               0.F,
0169                               &layer));
0170     LogDebug("LayerMeasurements") << "adding a missing hit.";
0171   } else
0172     LogDebug("LayerMeasurements") << "adding not measurement.";
0173
0174   return result;
0175 }
0176
0177 vector<TrajectoryMeasurementGroup> LayerMeasurements::groupedMeasurements(const DetLayer& layer,
0178                                                                           const TrajectoryStateOnSurface& startingState,
0179                                                                           const Propagator& prop,
0180                                                                           const MeasurementEstimator& est) const {
0181   vector<TrajectoryMeasurementGroup> result;
0182
0183   vector<DetGroup>&& groups = layer.groupedCompatibleDets(startingState, prop, est);
0184   result.reserve(groups.size());
0185
0186   tracking::TempMeasurements tmps;
0187   for (auto& grp : groups) {
0188     if (grp.empty())
0189       continue;
0190
0191     vector<TrajectoryMeasurement> tmpVec;
0192     for (auto const& det : grp) {
0193       auto const& detTS = det.trajectoryState();
0194
0195       if (not layerIsCompatibleWithTSOS(layer, detTS))
0196         continue;
0197
0198       MeasurementDetWithData mdet = detSystem_.idToDet(det.det()->geographicalId(), data_);
0199       if (mdet.isNull()) {
0200         throw MeasurementDetException("MeasurementDet not found");
0201       }
0202       if (mdet.measurements(detTS, est, tmps))
0203         for (std::size_t i = 0; i != tmps.size(); ++i)
0204           tmpVec.emplace_back(detTS, std::move(tmps.hits[i]), tmps.distances[i], &layer);
0205       tmps.clear();
0206     }
0207
0208     // sort the final result
0209     LogDebug("LayerMeasurements") << "Sorting " << tmpVec.size() << " measurements in this grp.";
0210     sort(tmpVec.begin(), tmpVec.end(), TrajMeasLessEstim());
0211     addInvalidMeas(tmpVec, grp, layer);
0212     result.emplace_back(std::move(tmpVec), std::move(grp));
0213   }
0214
0215   // if the result is empty check if the layer is compatible (for invalid measurement)
0216   if (result.empty()) {
0217     pair<bool, TrajectoryStateOnSurface> compat = layer.compatible(startingState, prop, est);
0218     if (compat.first and layerIsCompatibleWithTSOS(layer, compat.second)) {
0219       vector<TrajectoryMeasurement> tmVec;
0220       tmVec.emplace_back(compat.second,
0221                          std::make_shared<InvalidTrackingRecHitNoDet>(layer.surface(), TrackingRecHit::inactive),
0222                          0.F,
0223                          &layer);
0224       result.emplace_back(std::move(tmVec), DetGroup());
0225     }
0226   }
0227   return result;
0228 }