MeasurementDet/plugins/TkPhase2OTMeasurementDet.cc

0001 #include "TkPhase2OTMeasurementDet.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 "DataFormats/TrackerRecHit2D/interface/Phase2TrackerRecHit1D.h"
0007 #include "TrackingTools/DetLayers/interface/MeasurementEstimator.h"
0008 #include "TrackingTools/PatternTools/interface/TrajMeasLessEstim.h"
0009
0010 TkPhase2OTMeasurementDet::TkPhase2OTMeasurementDet(const GeomDet* gdet, Phase2OTMeasurementConditionSet& conditions)
0011     : MeasurementDet(gdet), index_(0), theDetConditions(&conditions) {
0012   if (dynamic_cast<const PixelGeomDetUnit*>(gdet) == nullptr) {
0013     throw MeasurementDetException(
0014         "TkPhase2OTMeasurementDet constructed with a GeomDet which is not a PixelGeomDetUnit");
0015   }
0016 }
0017
0018 bool TkPhase2OTMeasurementDet::measurements(const TrajectoryStateOnSurface& stateOnThisDet,
0019                                             const MeasurementEstimator& est,
0020                                             const MeasurementTrackerEvent& data,
0021                                             TempMeasurements& result) const {
0022   if (!isActive(data)) {
0023     result.add(theInactiveHit, 0.F);
0024     return true;
0025   }
0026
0027   if (recHits(stateOnThisDet, est, data, result.hits, result.distances))
0028     return true;
0029
0030   // create a TrajectoryMeasurement with an invalid RecHit and zero estimate
0031   bool inac = hasBadComponents(stateOnThisDet, data);
0032   result.add(inac ? theInactiveHit : theMissingHit, 0.F);
0033   return inac;
0034 }
0035
0036 bool TkPhase2OTMeasurementDet::recHits(const TrajectoryStateOnSurface& stateOnThisDet,
0037                                        const MeasurementEstimator& est,
0038                                        const MeasurementTrackerEvent& data,
0039                                        RecHitContainer& result,
0040                                        std::vector<float>& diffs) const {
0041   if UNLIKELY ((!isActive(data)) || isEmpty(data.phase2OTData()))
0042     return false;
0043
0044   auto oldSize = result.size();
0045
0046   const detset& detSet = data.phase2OTData().detSet(index());
0047   auto begin = &(data.phase2OTData().handle()->data().front());
0048   auto reject = [&](auto ci) -> bool {
0049     return (!data.phase2OTClustersToSkip().empty()) && data.phase2OTClustersToSkip()[ci - begin];
0050   };
0051
0052   /// use the usual 5 sigma cut from the Traj to identify the column....
0053   auto firstCluster = detSet.begin();
0054   auto lastCluster = detSet.end();
0055
0056   // do not use this as it does not account for APE...
0057   // auto xyLimits = est.maximalLocalDisplacement(stateOnThisDet,fastGeomDet().specificSurface());
0058   auto le = stateOnThisDet.localError().positionError();
0059   LocalError lape = static_cast<TrackerGeomDet const&>(fastGeomDet()).localAlignmentError();
0060   auto ye = le.yy();
0061   if (lape.valid()) {
0062     ye += lape.yy();
0063   }
0064   // 5 sigma to be on the safe side
0065   ye = 5.f * std::sqrt(ye);
0066   LocalVector maxD(0, ye, 0);
0067   // pixel topology is rectangular: x and y are independent
0068   auto ymin = specificGeomDet().specificTopology().measurementPosition(stateOnThisDet.localPosition() - maxD);
0069   auto ymax = specificGeomDet().specificTopology().measurementPosition(stateOnThisDet.localPosition() + maxD);
0070   int utraj = ymin.x();
0071   // do not apply for iteration not cutting on propagation
0072   if (est.maxSagitta() >= 0) {
0073     int colMin = ymin.y();
0074     int colMax = ymax.y();
0075     firstCluster = std::find_if(firstCluster, detSet.end(), [colMin](const Phase2TrackerCluster1D& hit) {
0076       return int(hit.column()) >= colMin;
0077     });
0078     lastCluster = std::find_if(
0079         firstCluster, detSet.end(), [colMax](const Phase2TrackerCluster1D& hit) { return int(hit.column()) > colMax; });
0080   }
0081
0082   while (firstCluster != lastCluster) {  // loop on each column
0083     auto const col = firstCluster->column();
0084     auto endCluster = std::find_if(
0085         firstCluster, detSet.end(), [col](const Phase2TrackerCluster1D& hit) { return hit.column() != col; });
0086     // find trajectory position in this column
0087     auto rightCluster = std::find_if(
0088         firstCluster, endCluster, [utraj](const Phase2TrackerCluster1D& hit) { return int(hit.firstStrip()) > utraj; });
0089     // search for compatible clusters...
0090     if (rightCluster != firstCluster) {
0091       // there are hits on the left of the utraj
0092       auto leftCluster = rightCluster;
0093       while (--leftCluster >= firstCluster) {
0094         if (reject(leftCluster))
0095           continue;
0096         Phase2TrackerCluster1DRef cluster = detSet.makeRefTo(data.phase2OTData().handle(), leftCluster);
0097         auto hit = buildRecHit(cluster, stateOnThisDet.localParameters());
0098         auto diffEst = est.estimate(stateOnThisDet, *hit);
0099         if (!diffEst.first)
0100           break;  // exit loop on first incompatible hit
0101         result.push_back(hit);
0102         diffs.push_back(diffEst.second);
0103       }
0104     }
0105     for (; rightCluster != endCluster; rightCluster++) {
0106       if (reject(rightCluster))
0107         continue;
0108       Phase2TrackerCluster1DRef cluster = detSet.makeRefTo(data.phase2OTData().handle(), rightCluster);
0109       auto hit = buildRecHit(cluster, stateOnThisDet.localParameters());
0110       auto diffEst = est.estimate(stateOnThisDet, *hit);
0111       if (!diffEst.first)
0112         break;  // exit loop on first incompatible hit
0113       result.push_back(hit);
0114       diffs.push_back(diffEst.second);
0115     }
0116     firstCluster = endCluster;
0117   }  // loop over columns
0118   return result.size() > oldSize;
0119 }
0120
0121 TrackingRecHit::RecHitPointer TkPhase2OTMeasurementDet::buildRecHit(const Phase2TrackerCluster1DRef& cluster,
0122                                                                     const LocalTrajectoryParameters& ltp) const {
0123   const PixelGeomDetUnit& gdu(specificGeomDet());
0124   auto&& params = cpe()->localParameters(*cluster, gdu);
0125
0126   return std::make_shared<Phase2TrackerRecHit1D>(params.first, params.second, fastGeomDet(), cluster);
0127 }
0128
0129 TkPhase2OTMeasurementDet::RecHitContainer TkPhase2OTMeasurementDet::recHits(const TrajectoryStateOnSurface& ts,
0130                                                                             const MeasurementTrackerEvent& data) const {
0131   RecHitContainer result;
0132   if (isEmpty(data.phase2OTData()))
0133     return result;
0134   if (!isActive(data))
0135     return result;
0136   const Phase2TrackerCluster1D* begin = nullptr;
0137   if (!data.phase2OTData().handle()->data().empty()) {
0138     begin = &(data.phase2OTData().handle()->data().front());
0139   }
0140   const detset& detSet = data.phase2OTData().detSet(index());
0141   result.reserve(detSet.size());
0142   for (const_iterator ci = detSet.begin(); ci != detSet.end(); ++ci) {
0143     if (ci < begin) {
0144       edm::LogError("IndexMisMatch") << "TkPhase2OTMeasurementDet cannot create hit because of index mismatch.";
0145       return result;
0146     }
0147     unsigned int index = ci - begin;
0148     if (!data.phase2OTClustersToSkip().empty() && index >= data.phase2OTClustersToSkip().size()) {
0149       edm::LogError("IndexMisMatch") << "TkPhase2OTMeasurementDet cannot create hit because of index mismatch. i.e "
0150                                      << index << " >= " << data.phase2OTClustersToSkip().size();
0151       return result;
0152     }
0153     if (data.phase2OTClustersToSkip().empty() or (not data.phase2OTClustersToSkip()[index])) {
0154       Phase2TrackerCluster1DRef cluster = detSet.makeRefTo(data.phase2OTData().handle(), ci);
0155       result.push_back(buildRecHit(cluster, ts.localParameters()));
0156     } else {
0157       LogDebug("TkPhase2OTMeasurementDet") << "skipping this cluster from last iteration on "
0158                                            << fastGeomDet().geographicalId().rawId() << " key: " << index;
0159     }
0160   }
0161   return result;
0162 }
0163
0164 //FIXME:just temporary solution for phase2!
0165 bool TkPhase2OTMeasurementDet::hasBadComponents(const TrajectoryStateOnSurface& tsos,
0166                                                 const MeasurementTrackerEvent& data) const {
0167   /*
0168     if (badRocPositions_.empty()) return false;
0169     LocalPoint lp = tsos.localPosition();
0170     LocalError le = tsos.localError().positionError();
0171     double dx = 3*std::sqrt(le.xx()) + theRocWidth, dy = 3*std::sqrt(le.yy()) + theRocHeight;
0172     for (std::vector<LocalPoint>::const_iterator it = badRocPositions_.begin(), ed = badRocPositions_.end(); it != ed; ++it) {
0173         if ( (std::abs(it->x() - lp.x()) < dx) &&
0174              (std::abs(it->y() - lp.y()) < dy) ) return true;
0175     }
0176 */
0177   return false;
0178 }