Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:28:12 

0001 #include "TkStripMeasurementDet.h"
0002 #include "TrackingTools/TransientTrackingRecHit/interface/InvalidTransientRecHit.h"
0003 #include "Geometry/CommonTopologies/interface/StripTopology.h"
0004 #include "TrackingTools/MeasurementDet/interface/MeasurementDetException.h"
0005 #include "TrackingTools/PatternTools/interface/TrajectoryMeasurement.h"
0006 #include "RecoTracker/TransientTrackingRecHit/interface/TSiStripRecHit2DLocalPos.h"
0007 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h"
0008 #include "TrackingTools/DetLayers/interface/MeasurementEstimator.h"
0009 #include "TrackingTools/PatternTools/interface/TrajMeasLessEstim.h"
0010 
0011 #include <typeinfo>
0012 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0013 #include "DataFormats/SiStripDetId/interface/SiStripDetId.h"
0014 #include <cassert>
0015 
0016 TkStripMeasurementDet::TkStripMeasurementDet(const GeomDet* gdet, StMeasurementConditionSet& conditions)
0017     : MeasurementDet(gdet), index_(-1), theDetConditions(&conditions) {
0018   if (dynamic_cast<const StripGeomDetUnit*>(gdet) == nullptr) {
0019     throw MeasurementDetException("TkStripMeasurementDet constructed with a GeomDet which is not a StripGeomDetUnit");
0020   }
0021 }
0022 
0023 // fast check if the det contains any useful cluster
0024 bool TkStripMeasurementDet::empty(const MeasurementTrackerEvent& data) const {
0025   if UNLIKELY ((!isActive(data)) || isEmpty(data.stripData()))
0026     return true;
0027 
0028   const detset& detSet = data.stripData().detSet(index());
0029   for (auto ci = detSet.begin(); ci != detSet.end(); ++ci) {
0030     if (isMasked(*ci))
0031       continue;
0032     SiStripClusterRef cluster = detSet.makeRefTo(data.stripData().handle(), ci);
0033     if (accept(cluster, data.stripClustersToSkip()))
0034       return false;
0035   }
0036   return true;
0037 }
0038 
0039 TkStripMeasurementDet::RecHitContainer TkStripMeasurementDet::recHits(const TrajectoryStateOnSurface& ts,
0040                                                                       const MeasurementTrackerEvent& data) const {
0041   RecHitContainer result;
0042   if UNLIKELY ((!isActive(data)) || isEmpty(data.stripData()))
0043     return result;
0044   const detset& detSet = data.stripData().detSet(index());
0045   result.reserve(detSet.size());
0046   for (new_const_iterator ci = detSet.begin(); ci != detSet.end(); ++ci) {
0047     if (isMasked(*ci))
0048       continue;
0049     // for ( ClusterIterator ci=theClusterRange.first; ci != theClusterRange.second; ci++) {
0050     SiStripClusterRef cluster = detSet.makeRefTo(data.stripData().handle(), ci);
0051     if (accept(cluster, data.stripClustersToSkip()))
0052       result.push_back(buildRecHit(cluster, ts));
0053     else
0054       LogDebug("TkStripMeasurementDet") << "skipping this str from last iteration on" << rawId()
0055                                         << " key: " << cluster.key();
0056   }
0057   return result;
0058 }
0059 
0060 // FIXME need to be merged with simpleRecHits
0061 bool TkStripMeasurementDet::recHits(SimpleHitContainer& result,
0062                                     const TrajectoryStateOnSurface& stateOnThisDet,
0063                                     const MeasurementEstimator& est,
0064                                     const MeasurementTrackerEvent& data) const {
0065   if UNLIKELY ((!isActive(data)) || isEmpty(data.stripData()))
0066     return false;
0067   auto oldSize = result.size();
0068 
0069   int utraj = specificGeomDet().specificTopology().measurementPosition(stateOnThisDet.localPosition()).x();
0070   const detset& detSet = data.stripData().detSet(index());
0071   auto const& cpepar = cpe()->getAlgoParam(specificGeomDet(), stateOnThisDet.localParameters());
0072 
0073   auto rightCluster = std::find_if(
0074       detSet.begin(), detSet.end(), [utraj](const SiStripCluster& hit) { return hit.firstStrip() > utraj; });
0075 
0076   std::vector<SiStripRecHit2D> tmp;
0077   if (rightCluster != detSet.begin()) {
0078     // there are hits on the left of the utraj
0079     auto leftCluster = rightCluster;
0080     while (--leftCluster >= detSet.begin()) {
0081       SiStripClusterRef clusterref = detSet.makeRefTo(data.stripData().handle(), leftCluster);
0082       bool isCompatible = filteredRecHits(clusterref, cpepar, stateOnThisDet, est, data.stripClustersToSkip(), tmp);
0083       if (!isCompatible)
0084         break;  // exit loop on first incompatible hit
0085       for (auto&& h : tmp)
0086         result.push_back(new SiStripRecHit2D(h));
0087       tmp.clear();
0088     }
0089   }
0090   for (; rightCluster != detSet.end(); rightCluster++) {
0091     SiStripClusterRef clusterref = detSet.makeRefTo(data.stripData().handle(), rightCluster);
0092     bool isCompatible = filteredRecHits(clusterref, cpepar, stateOnThisDet, est, data.stripClustersToSkip(), tmp);
0093     if (!isCompatible)
0094       break;  // exit loop on first incompatible hit
0095     for (auto&& h : tmp)
0096       result.push_back(new SiStripRecHit2D(h));
0097     tmp.clear();
0098   }
0099 
0100   return result.size() > oldSize;
0101 }
0102 
0103 bool TkStripMeasurementDet::simpleRecHits(const TrajectoryStateOnSurface& stateOnThisDet,
0104                                           const MeasurementEstimator& est,
0105                                           const MeasurementTrackerEvent& data,
0106                                           std::vector<SiStripRecHit2D>& result) const {
0107   if UNLIKELY ((!isActive(data)) || isEmpty(data.stripData()))
0108     return false;
0109 
0110   auto oldSize = result.size();
0111 
0112   int utraj = specificGeomDet().specificTopology().measurementPosition(stateOnThisDet.localPosition()).x();
0113   const detset& detSet = data.stripData().detSet(index());
0114   auto const& cpepar = cpe()->getAlgoParam(specificGeomDet(), stateOnThisDet.localParameters());
0115 
0116   auto rightCluster = std::find_if(
0117       detSet.begin(), detSet.end(), [utraj](const SiStripCluster& hit) { return hit.firstStrip() > utraj; });
0118 
0119   if (rightCluster != detSet.begin()) {
0120     // there are hits on the left of the utraj
0121     auto leftCluster = rightCluster;
0122     while (--leftCluster >= detSet.begin()) {
0123       SiStripClusterRef clusterref = detSet.makeRefTo(data.stripData().handle(), leftCluster);
0124       bool isCompatible = filteredRecHits(clusterref, cpepar, stateOnThisDet, est, data.stripClustersToSkip(), result);
0125       if (!isCompatible)
0126         break;  // exit loop on first incompatible hit
0127     }
0128   }
0129   for (; rightCluster != detSet.end(); rightCluster++) {
0130     SiStripClusterRef clusterref = detSet.makeRefTo(data.stripData().handle(), rightCluster);
0131     bool isCompatible = filteredRecHits(clusterref, cpepar, stateOnThisDet, est, data.stripClustersToSkip(), result);
0132     if (!isCompatible)
0133       break;  // exit loop on first incompatible hit
0134   }
0135 
0136   return result.size() > oldSize;
0137 }
0138 
0139 bool TkStripMeasurementDet::recHits(const TrajectoryStateOnSurface& stateOnThisDet,
0140                                     const MeasurementEstimator& est,
0141                                     const MeasurementTrackerEvent& data,
0142                                     RecHitContainer& result,
0143                                     std::vector<float>& diffs) const {
0144   if UNLIKELY ((!isActive(data)) || isEmpty(data.stripData()))
0145     return false;
0146 
0147   auto oldSize = result.size();
0148   auto const& cpepar = cpe()->getAlgoParam(specificGeomDet(), stateOnThisDet.localParameters());
0149 
0150   int utraj = specificGeomDet().specificTopology().measurementPosition(stateOnThisDet.localPosition()).x();
0151 
0152   const detset& detSet = data.stripData().detSet(index());
0153   auto rightCluster = std::find_if(
0154       detSet.begin(), detSet.end(), [utraj](const SiStripCluster& hit) { return hit.firstStrip() > utraj; });
0155 
0156   if (rightCluster != detSet.begin()) {
0157     // there are hits on the left of the utraj
0158     auto leftCluster = rightCluster;
0159     while (--leftCluster >= detSet.begin()) {
0160       SiStripClusterRef clusterref = detSet.makeRefTo(data.stripData().handle(), leftCluster);
0161       bool isCompatible =
0162           filteredRecHits(clusterref, cpepar, stateOnThisDet, est, data.stripClustersToSkip(), result, diffs);
0163       if (!isCompatible)
0164         break;  // exit loop on first incompatible hit
0165     }
0166   }
0167   for (; rightCluster != detSet.end(); rightCluster++) {
0168     SiStripClusterRef clusterref = detSet.makeRefTo(data.stripData().handle(), rightCluster);
0169     bool isCompatible =
0170         filteredRecHits(clusterref, cpepar, stateOnThisDet, est, data.stripClustersToSkip(), result, diffs);
0171     if (!isCompatible)
0172       break;  // exit loop on first incompatible hit
0173   }
0174 
0175   return result.size() > oldSize;
0176 }
0177 
0178 bool TkStripMeasurementDet::measurements(const TrajectoryStateOnSurface& stateOnThisDet,
0179                                          const MeasurementEstimator& est,
0180                                          const MeasurementTrackerEvent& data,
0181                                          TempMeasurements& result) const {
0182   if (!isActive(data)) {
0183     LogDebug("TkStripMeasurementDet") << " found an inactive module " << rawId();
0184     result.add(theInactiveHit, 0.F);
0185     return true;
0186   }
0187 
0188   if (!isEmpty(data.stripData())) {
0189     LogDebug("TkStripMeasurementDet") << " found hit on this module " << rawId();
0190     RecHitContainer rechits;
0191     std::vector<float> diffs;
0192     if (recHits(stateOnThisDet, est, data, result.hits, result.distances))
0193       return true;
0194   }
0195 
0196   // create a TrajectoryMeasurement with an invalid RecHit and zero estimate
0197 
0198   if (!stateOnThisDet.hasError()) {
0199     result.add(theMissingHit, 0.F);
0200     return false;
0201   }
0202 
0203   float utraj = specificGeomDet().specificTopology().measurementPosition(stateOnThisDet.localPosition()).x();
0204   float uerr = sqrt(specificGeomDet()
0205                         .specificTopology()
0206                         .measurementError(stateOnThisDet.localPosition(), stateOnThisDet.localError().positionError())
0207                         .uu());
0208   if (testStrips(utraj, uerr)) {
0209     //LogDebug("TkStripMeasurementDet") << " DetID " << id_ << " empty after search, but active ";
0210     result.add(theMissingHit, 0.F);
0211     return false;
0212   }
0213 
0214   //LogDebug("TkStripMeasurementDet") << " DetID " << id_ << " empty after search, and inactive ";
0215   result.add(theInactiveHit, 0.F);
0216   return true;
0217 }
0218 
0219 void TkStripMeasurementDet::simpleRecHits(const TrajectoryStateOnSurface& ts,
0220                                           const MeasurementTrackerEvent& data,
0221                                           std::vector<SiStripRecHit2D>& result) const {
0222   if (isEmpty(data.stripData()) || !isActive(data))
0223     return;
0224 
0225   const detset& detSet = data.stripData().detSet(index());
0226   unInitDynArray(AClusters::value_type, detSet.size(), clusters);
0227   assert(clusters.empty());
0228   for (auto const& ci : detSet) {
0229     if (isMasked(ci))
0230       continue;
0231     if (accept(detSet.makeKeyOf(&ci), data.stripClustersToSkip()))
0232       clusters.push_back(&ci);
0233     else
0234       LogDebug("TkStripMeasurementDet") << "skipping this str from last iteration on" << rawId()
0235                                         << " key: " << detSet.makeKeyOf(&ci);
0236   }
0237   if (!clusters.empty())
0238     buildSimpleRecHits(clusters, data, detSet, ts, result);
0239 }
0240 
0241 std::tuple<TkStripRecHitIter, TkStripRecHitIter> TkStripMeasurementDet::hitRange(
0242     const TrajectoryStateOnSurface& ts, const MeasurementTrackerEvent& data) const {
0243   if (isEmpty(data.stripData()) || !isActive(data))
0244     return std::tuple<TkStripRecHitIter, TkStripRecHitIter>();
0245   const detset& detSet = data.stripData().detSet(index());
0246   return std::make_tuple(TkStripRecHitIter(detSet.begin(), detSet.end(), *this, ts, data),
0247                          TkStripRecHitIter(detSet.end(), detSet.end(), *this, ts, data));
0248 }
0249 
0250 void TkStripMeasurementDet::advance(TkStripRecHitIter& hi) const {
0251   while (!hi.empty()) {
0252     auto ci = hi.clusterI;
0253     auto const& data = *hi.data;
0254     if (isMasked(*ci))
0255       continue;
0256     SiStripClusterRef cluster = edmNew::makeRefTo(data.stripData().handle(), ci);
0257     if (accept(cluster, data.stripClustersToSkip()))
0258       return;
0259     ++hi.clusterI;
0260   }
0261 }
0262 
0263 SiStripRecHit2D TkStripMeasurementDet::hit(TkStripRecHitIter const& hi) const {
0264   const GeomDetUnit& gdu(specificGeomDet());
0265   auto ci = hi.clusterI;
0266   auto const& data = *hi.data;
0267   auto const& ltp = *hi.tsos;
0268 
0269   SiStripClusterRef cluster = edmNew::makeRefTo(data.stripData().handle(), ci);
0270   LocalValues lv = cpe()->localParameters(*cluster, gdu, ltp);
0271   return SiStripRecHit2D(lv.first, lv.second, gdu, cluster);
0272 }
0273 
0274 bool TkStripMeasurementDet::testStrips(float utraj, float uerr) const {
0275   int16_t start = (int16_t)std::max<float>(utraj - 3.f * uerr, 0);
0276   int16_t end = (int16_t)std::min<float>(utraj + 3.f * uerr, totalStrips());
0277 
0278   if (start >= end) {  // which means either end <=0 or start >= totalStrips_
0279     /* LogDebug("TkStripMeasurementDet") << "Testing module " << id_ <<","<<
0280             " U = " << utraj << " +/- " << uerr << 
0281             "; Range [" << (utraj - 3*uerr) << ", " << (utraj + 3*uerr) << "] " << 
0282             ": YOU'RE COMPLETELY OFF THE MODULE."; */
0283     //return false;
0284     return true;  // Wolfgang thinks this way is better
0285                   // and solves some problems with grouped ckf
0286   }
0287 
0288   typedef std::vector<BadStripBlock>::const_iterator BSBIT;
0289 
0290   int16_t bad = 0, largestBadBlock = 0;
0291   for (BSBIT bsbc = badStripBlocks().begin(), bsbe = badStripBlocks().end(); bsbc != bsbe; ++bsbc) {
0292     if (bsbc->last < start)
0293       continue;
0294     if (bsbc->first > end)
0295       break;
0296     int16_t thisBad = std::min(bsbc->last, end) - std::max(bsbc->first, start);
0297     if (thisBad > largestBadBlock)
0298       largestBadBlock = thisBad;
0299     bad += thisBad;
0300   }
0301 
0302   bool ok = (bad < (end - start)) && (uint16_t(bad) <= badStripCuts().maxBad) &&
0303             (uint16_t(largestBadBlock) <= badStripCuts().maxConsecutiveBad);
0304 
0305   //    if (bad) {
0306   //       edm::LogWarning("TkStripMeasurementDet") << "Testing module " << id_ <<" (subdet: "<< SiStripDetId(id_).subdetId() << ")" <<
0307   //            " U = " << utraj << " +/- " << uerr <<
0308   //            "; Range [" << (utraj - 3*uerr) << ", " << (utraj + 3*uerr) << "] " <<
0309   //            "= [" << start << "," << end << "]" <<
0310   //            " total strips:" << (end-start) << ", good:" << (end-start-bad) << ", bad:" << bad << ", largestBadBlock: " << largestBadBlock <<
0311   //            ". " << (ok ? "OK" : "NO");
0312   //    }
0313   return ok;
0314 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!