MeasurementDet/plugins/TkPixelMeasurementDet.cc

0001 #include "TkPixelMeasurementDet.h"
0002 #include "TrackingTools/TransientTrackingRecHit/interface/InvalidTransientRecHit.h"
0003 #include "Geometry/CommonTopologies/interface/PixelTopology.h"
0004 #include "TrackingTools/MeasurementDet/interface/MeasurementDetException.h"
0005 #include "TrackingTools/PatternTools/interface/TrajectoryMeasurement.h"
0006 #include "RecoTracker/TransientTrackingRecHit/interface/TSiPixelRecHit.h"
0007 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
0008 #include "TrackingTools/DetLayers/interface/MeasurementEstimator.h"
0009 #include "TrackingTools/PatternTools/interface/TrajMeasLessEstim.h"
0010
0011 namespace {
0012   // in cms units are in cm
0013   constexpr float theRocWidth = 0.81 / 2;
0014   constexpr float theRocHeight = 0.81 / 2;
0015 }  // namespace
0016
0017 TkPixelMeasurementDet::TkPixelMeasurementDet(const GeomDet* gdet, PxMeasurementConditionSet& conditions)
0018     : MeasurementDet(gdet), index_(0), theDetConditions(&conditions) {
0019   if (dynamic_cast<const PixelGeomDetUnit*>(gdet) == nullptr) {
0020     throw MeasurementDetException("TkPixelMeasurementDet constructed with a GeomDet which is not a PixelGeomDetUnit");
0021   }
0022 }
0023
0024 bool TkPixelMeasurementDet::measurements(const TrajectoryStateOnSurface& stateOnThisDet,
0025                                          const MeasurementEstimator& est,
0026                                          const MeasurementTrackerEvent& data,
0027                                          TempMeasurements& result) const {
0028   if (!isActive(data)) {
0029     result.add(theInactiveHit, 0.F);
0030     return true;
0031   }
0032
0033   auto xl = 100.f;
0034   auto yl = 100.f;
0035   // do not apply for iteration not cutting on propagation
0036   if (est.maxSagitta() >= 0) {
0037     // do not use this as it does not account for APE...
0038     // auto xyLimits = est.maximalLocalDisplacement(stateOnThisDet,fastGeomDet().specificSurface());
0039     auto le = stateOnThisDet.localError().positionError();
0040     LocalError lape = static_cast<TrackerGeomDet const&>(fastGeomDet()).localAlignmentError();
0041     xl = le.xx();
0042     yl = le.yy();
0043     if (lape.valid()) {
0044       xl += lape.xx();
0045       yl += lape.yy();
0046     }
0047     // 5 sigma to be on the safe side
0048     xl = 5.f * std::sqrt(xl);
0049     yl = 5.f * std::sqrt(yl);
0050   }
0051
0052   auto oldSize = result.size();
0053   MeasurementDet::RecHitContainer&& allHits = compHits(stateOnThisDet, data, xl, yl);
0054   for (auto&& hit : allHits) {
0055     std::pair<bool, double> diffEst = est.estimate(stateOnThisDet, *hit);
0056     if (diffEst.first)
0057       result.add(std::move(hit), diffEst.second);
0058   }
0059
0060   if (result.size() > oldSize)
0061     return true;
0062
0063   // create a TrajectoryMeasurement with an invalid RecHit and zero estimate
0064   bool inac = hasBadComponents(stateOnThisDet, data);
0065   result.add(inac ? theInactiveHit : theMissingHit, 0.F);
0066   return inac;
0067 }
0068
0069 TrackingRecHit::RecHitPointer TkPixelMeasurementDet::buildRecHit(const SiPixelClusterRef& cluster,
0070                                                                  const LocalTrajectoryParameters& ltp) const {
0071   const GeomDetUnit& gdu(specificGeomDet());
0072
0073   auto&& params = cpe()->getParameters(*cluster, gdu, ltp);
0074   return std::make_shared<SiPixelRecHit>(
0075       std::get<0>(params), std::get<1>(params), std::get<2>(params), fastGeomDet(), cluster);
0076 }
0077
0078 TkPixelMeasurementDet::RecHitContainer TkPixelMeasurementDet::recHits(const TrajectoryStateOnSurface& ts,
0079                                                                       const MeasurementTrackerEvent& data) const {
0080   float xl = 100.f;  // larger than any detector
0081   float yl = 100.f;
0082   return compHits(ts, data, xl, yl);
0083 }
0084
0085 TkPixelMeasurementDet::RecHitContainer TkPixelMeasurementDet::compHits(const TrajectoryStateOnSurface& ts,
0086                                                                        const MeasurementTrackerEvent& data,
0087                                                                        float xl,
0088                                                                        float yl) const {
0089   RecHitContainer result;
0090   if (isEmpty(data.pixelData()) == true)
0091     return result;
0092   if (isActive(data) == false)
0093     return result;
0094   const SiPixelCluster* begin = nullptr;
0095   if (!data.pixelData().handle()->data().empty()) {
0096     begin = &(data.pixelData().handle()->data().front());
0097   }
0098   const detset& detSet = data.pixelData().detSet(index());
0099   result.reserve(detSet.size());
0100
0101   // pixel topology is rectangular, all positions are independent
0102   LocalVector maxD(xl, yl, 0);
0103   auto PMinus = specificGeomDet().specificTopology().measurementPosition(ts.localPosition() - maxD);
0104   auto PPlus = specificGeomDet().specificTopology().measurementPosition(ts.localPosition() + maxD);
0105
0106   int xminus = PMinus.x();
0107   int yminus = PMinus.y();
0108   int xplus = PPlus.x() + 0.5f;
0109   int yplus = PPlus.y() + 0.5f;
0110
0111   // rechits are sorted in x...
0112   auto rightCluster = std::find_if(
0113       detSet.begin(), detSet.end(), [xplus](const SiPixelCluster& cl) { return cl.minPixelRow() > xplus; });
0114
0115   // std::cout << "px xlim " << xl << ' ' << xminus << '/' << xplus << ' ' << rightCluster-detSet.begin() << ',' << detSet.end()-rightCluster << std::endl;
0116
0117   // consider only compatible clusters
0118   for (auto ci = detSet.begin(); ci != rightCluster; ++ci) {
0119     if (ci < begin) {
0120       edm::LogError("IndexMisMatch") << "TkPixelMeasurementDet cannot create hit because of index mismatch.";
0121       return result;
0122     }
0123     unsigned int index = ci - begin;
0124     if (!data.pixelClustersToSkip().empty() && index >= data.pixelClustersToSkip().size()) {
0125       edm::LogError("IndexMisMatch") << "TkPixelMeasurementDet cannot create hit because of index mismatch. i.e "
0126                                      << index << " >= " << data.pixelClustersToSkip().size();
0127       return result;
0128     }
0129
0130     if (ci->maxPixelRow() < xminus)
0131       continue;
0132     // also check compatibility in y... (does not add much)
0133     if (ci->minPixelCol() > yplus)
0134       continue;
0135     if (ci->maxPixelCol() < yminus)
0136       continue;
0137
0138     if (data.pixelClustersToSkip().empty() or (not data.pixelClustersToSkip()[index])) {
0139       SiPixelClusterRef cluster = detSet.makeRefTo(data.pixelData().handle(), ci);
0140       result.push_back(buildRecHit(cluster, ts.localParameters()));
0141     } else {
0142       LogDebug("TkPixelMeasurementDet") << "skipping this cluster from last iteration on "
0143                                         << fastGeomDet().geographicalId().rawId() << " key: " << index;
0144     }
0145   }
0146   return result;
0147 }
0148
0149 bool TkPixelMeasurementDet::hasBadComponents(const TrajectoryStateOnSurface& tsos,
0150                                              const MeasurementTrackerEvent& data) const {
0151   auto badFEDChannelPositions = getBadFEDChannelPositions(data);
0152   if (badRocPositions_.empty() && badFEDChannelPositions == nullptr)
0153     return false;
0154
0155   auto lp = tsos.localPosition();
0156   auto le = tsos.localError().positionError();
0157   for (auto const& broc : badRocPositions_) {
0158     auto dx = std::abs(broc.x() - lp.x()) - theRocWidth;
0159     auto dy = std::abs(broc.y() - lp.y()) - theRocHeight;
0160     if ((dx <= 0.f) & (dy <= 0.f))
0161       return true;
0162     if ((dx * dx < 9.f * le.xx()) && (dy * dy < 9.f * le.yy()))
0163       return true;
0164   }
0165
0166   if (badFEDChannelPositions == nullptr)
0167     return false;
0168   float dx = 3.f * std::sqrt(le.xx()) + theRocWidth, dy = 3.f * std::sqrt(le.yy()) + theRocHeight;
0169   for (auto const& p : *badFEDChannelPositions) {
0170     if (lp.x() > (p.first.x() - dx) && lp.x() < (p.second.x() + dx) && lp.y() > (p.first.y() - dy) &&
0171         lp.y() < (p.second.y() + dy)) {
0172       return true;
0173     }
0174   }
0175
0176   return false;
0177 }