PixelLowPtUtilities/src/StripSubClusterShapeTrajectoryFilter.cc

0001 #include "RecoTracker/PixelLowPtUtilities/interface/StripSubClusterShapeTrajectoryFilter.h"
0002
0003 #include <map>
0004 #include <set>
0005 #include <algorithm>
0006
0007 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0008 #include "DataFormats/TrackerRecHit2D/interface/ProjectedSiStripRecHit2D.h"
0009 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
0010 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
0011 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h"
0012 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
0013 #include "DataFormats/SiStripCommon/interface/ConstantsForHardwareSystems.h"
0014 #include "FWCore/Framework/interface/ESHandle.h"
0015 #include "FWCore/Framework/interface/EventSetup.h"
0016 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0017 #include "FWCore/ServiceRegistry/interface/Service.h"
0018 #include "Geometry/CommonDetUnit/interface/GeomDetType.h"
0019 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0020 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
0021 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0022 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0023 #include "MagneticField/Engine/interface/MagneticField.h"
0024 #include "RecoTracker/PixelLowPtUtilities/interface/ClusterShapeHitFilter.h"
0025 #include "RecoTracker/Record/interface/CkfComponentsRecord.h"
0026 #include "RecoTracker/MeasurementDet/interface/MeasurementTrackerEvent.h"
0027 #include "TrackingTools/PatternTools/interface/TempTrajectory.h"
0028 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0029 #include "FWCore/Framework/interface/ConsumesCollector.h"
0030 #include "RecoTracker/TkSeedGenerator/interface/FastHelix.h"
0031 #include "RecoTracker/TkSeedingLayers/interface/SeedingHitSet.h"
0032 #include "DataFormats/SiStripDetId/interface/SiStripDetId.h"
0033
0034 #ifdef StripSubClusterShapeFilterBase_COUNTERS
0035 #define INC_COUNTER(X) X++;
0036 #else
0037 #define INC_COUNTER(X)
0038 #endif
0039 namespace {
0040   class SlidingPeakFinder {
0041   public:
0042     SlidingPeakFinder(unsigned int size) : size_(size), half_((size + 1) / 2) {}
0043
0044     template <typename Test>
0045     bool apply(const uint8_t *x,
0046                const uint8_t *begin,
0047                const uint8_t *end,
0048                const Test &test,
0049                bool verbose = false,
0050                int firststrip = 0) {
0051       const uint8_t *ileft = (x != begin) ? std::min_element(x - 1, x + half_) : begin - 1;
0052       const uint8_t *iright = ((x + size_) < end) ? std::min_element(x + half_, std::min(x + size_ + 1, end)) : end;
0053       uint8_t left = (ileft < begin ? 0 : *ileft);
0054       uint8_t right = (iright >= end ? 0 : *iright);
0055       uint8_t center = *std::max_element(x, std::min(x + size_, end));
0056       uint8_t maxmin = std::max(left, right);
0057       if (maxmin < center) {
0058         bool ret = test(center, maxmin);
0059         if (ret) {
0060           ret = test(ileft, iright, begin, end);
0061         }
0062         return ret;
0063       } else {
0064         return false;
0065       }
0066     }
0067
0068     template <typename V, typename Test>
0069     bool apply(const V &ampls, const Test &test, bool verbose = false, int firststrip = 0) {
0070       const uint8_t *begin = &*ampls.begin();
0071       const uint8_t *end = &*ampls.end();
0072       for (const uint8_t *x = begin; x < end - (half_ - 1); ++x) {
0073         if (apply(x, begin, end, test, verbose, firststrip)) {
0074           return true;
0075         }
0076       }
0077       return false;
0078     }
0079
0080   private:
0081     unsigned int size_, half_;
0082   };
0083
0084   struct PeakFinderTest {
0085     PeakFinderTest(float mip,
0086                    uint32_t detid,
0087                    uint32_t firstStrip,
0088                    const SiStripNoises *theNoise,
0089                    float seedCutMIPs,
0090                    float seedCutSN,
0091                    float subclusterCutMIPs,
0092                    float subclusterCutSN)
0093         : mip_(mip),
0094           detid_(detid),
0095           firstStrip_(firstStrip),
0096           noiseObj_(theNoise),
0097           noises_(theNoise->getRange(detid)),
0098           subclusterCutMIPs_(subclusterCutMIPs),
0099           sumCut_(subclusterCutMIPs_ * mip_),
0100           subclusterCutSN2_(subclusterCutSN * subclusterCutSN) {
0101       cut_ = std::min<float>(seedCutMIPs * mip, seedCutSN * noiseObj_->getNoise(firstStrip + 1, noises_));
0102     }
0103
0104     bool operator()(uint8_t max, uint8_t min) const { return max - min > cut_; }
0105     bool operator()(const uint8_t *left, const uint8_t *right, const uint8_t *begin, const uint8_t *end) const {
0106       int yleft = (left < begin ? 0 : *left);
0107       int yright = (right >= end ? 0 : *right);
0108       float sum = 0.0;
0109       int maxval = 0;
0110       float noise = 0;
0111       for (const uint8_t *x = left + 1; x < right; ++x) {
0112         int baseline = (yleft * int(right - x) + yright * int(x - left)) / int(right - left);
0113         sum += int(*x) - baseline;
0114         noise += std::pow(noiseObj_->getNoise(firstStrip_ + int(x - begin), noises_), 2);
0115         maxval = std::max(maxval, int(*x) - baseline);
0116       }
0117       if (sum > sumCut_ && sum * sum > noise * subclusterCutSN2_)
0118         return true;
0119       return false;
0120     }
0121
0122   private:
0123     float mip_;
0124     unsigned int detid_;
0125     int firstStrip_;
0126     const SiStripNoises *noiseObj_;
0127     SiStripNoises::Range noises_;
0128     uint8_t cut_;
0129     float subclusterCutMIPs_, sumCut_, subclusterCutSN2_;
0130   };
0131
0132 }  // namespace
0133
0134 StripSubClusterShapeFilterBase::StripSubClusterShapeFilterBase(const edm::ParameterSet &iCfg,
0135                                                                edm::ConsumesCollector &iC)
0136     : topoToken_(iC.esConsumes<TrackerTopology, TrackerTopologyRcd>()),
0137       csfToken_(
0138           iC.esConsumes<ClusterShapeHitFilter, CkfComponentsRecord>(edm::ESInputTag("", "ClusterShapeHitFilter"))),
0139       geomToken_(iC.esConsumes<TrackerGeometry, TrackerDigiGeometryRecord>()),
0140       stripNoiseToken_(iC.esConsumes<SiStripNoises, SiStripNoisesRcd>()),
0141       label_(iCfg.getUntrackedParameter<std::string>("label", "")),
0142       maxNSat_(iCfg.getParameter<uint32_t>("maxNSat")),
0143       trimMaxADC_(iCfg.getParameter<double>("trimMaxADC")),
0144       trimMaxFracTotal_(iCfg.getParameter<double>("trimMaxFracTotal")),
0145       trimMaxFracNeigh_(iCfg.getParameter<double>("trimMaxFracNeigh")),
0146       maxTrimmedSizeDiffPos_(iCfg.getParameter<double>("maxTrimmedSizeDiffPos")),
0147       maxTrimmedSizeDiffNeg_(iCfg.getParameter<double>("maxTrimmedSizeDiffNeg")),
0148       subclusterWindow_(iCfg.getParameter<double>("subclusterWindow")),
0149       seedCutMIPs_(iCfg.getParameter<double>("seedCutMIPs")),
0150       seedCutSN_(iCfg.getParameter<double>("seedCutSN")),
0151       subclusterCutMIPs_(iCfg.getParameter<double>("subclusterCutMIPs")),
0152       subclusterCutSN_(iCfg.getParameter<double>("subclusterCutSN"))
0153 #ifdef StripSubClusterShapeFilterBase_COUNTERS
0154       ,
0155       called_(0),
0156       saturated_(0),
0157       test_(0),
0158       passTrim_(0),
0159       failTooLarge_(0),
0160       passSC_(0),
0161       failTooNarrow_(0)
0162 #endif
0163 {
0164   const edm::ParameterSet &iLM = iCfg.getParameter<edm::ParameterSet>("layerMask");
0165   if (not iLM.empty()) {
0166     const char *ndets[4] = {"TIB", "TID", "TOB", "TEC"};
0167     const int idets[4] = {3, 4, 5, 6};
0168     for (unsigned int i = 0; i < 4; ++i) {
0169       if (iLM.existsAs<bool>(ndets[i])) {
0170         std::fill(layerMask_[idets[i]].begin(), layerMask_[idets[i]].end(), iLM.getParameter<bool>(ndets[i]));
0171       } else {
0172         layerMask_[idets[i]][0] = 2;
0173         std::fill(layerMask_[idets[i]].begin() + 1, layerMask_[idets[i]].end(), 0);
0174         for (uint32_t lay : iLM.getParameter<std::vector<uint32_t>>(ndets[i])) {
0175           layerMask_[idets[i]][lay] = 1;
0176         }
0177       }
0178     }
0179   } else {
0180     for (auto &arr : layerMask_) {
0181       std::fill(arr.begin(), arr.end(), 1);
0182     }
0183   }
0184 }
0185
0186 StripSubClusterShapeFilterBase::~StripSubClusterShapeFilterBase() {
0187 #if 0
0188   std::cout << "StripSubClusterShapeFilterBase " << label_ << ": called        " << called_ << std::endl;
0189   std::cout << "StripSubClusterShapeFilterBase " << label_ << ": saturated     " << saturated_ << std::endl;
0190   std::cout << "StripSubClusterShapeFilterBase " << label_ << ": test          " << test_ << std::endl;
0191   std::cout << "StripSubClusterShapeFilterBase " << label_ << ": failTooNarrow " << failTooNarrow_ << std::endl;
0192   std::cout << "StripSubClusterShapeFilterBase " << label_ << ": passTrim      " << passTrim_ << std::endl;
0193   std::cout << "StripSubClusterShapeFilterBase " << label_ << ": passSC        " << passSC_ << std::endl;
0194   std::cout << "StripSubClusterShapeFilterBase " << label_ << ": failTooLarge  " << failTooLarge_ << std::endl;
0195 #endif
0196 }
0197
0198 void StripSubClusterShapeFilterBase::fillPSetDescription(edm::ParameterSetDescription &iDesc) {
0199   iDesc.addUntracked<std::string>("label", "");
0200   iDesc.add<uint32_t>("maxNSat", 3);
0201   iDesc.add<double>("trimMaxADC", 30.);
0202   iDesc.add<double>("trimMaxFracTotal", .15);
0203   iDesc.add<double>("trimMaxFracNeigh", .25);
0204   iDesc.add<double>("maxTrimmedSizeDiffPos", .7);
0205   iDesc.add<double>("maxTrimmedSizeDiffNeg", 1.);
0206   iDesc.add<double>("subclusterWindow", .7);
0207   iDesc.add<double>("seedCutMIPs", .35);
0208   iDesc.add<double>("seedCutSN", 7.);
0209   iDesc.add<double>("subclusterCutMIPs", .45);
0210   iDesc.add<double>("subclusterCutSN", 12.);
0211
0212   edm::ParameterSetDescription psdLM;
0213   psdLM.setAllowAnything();
0214   iDesc.add<edm::ParameterSetDescription>("layerMask", psdLM);
0215 }
0216
0217 bool StripSubClusterShapeFilterBase::testLastHit(const TrackingRecHit *hit,
0218                                                  const TrajectoryStateOnSurface &tsos,
0219                                                  bool mustProject) const {
0220   return testLastHit(hit, tsos.globalPosition(), tsos.globalDirection(), mustProject);
0221 }
0222 bool StripSubClusterShapeFilterBase::testLastHit(const TrackingRecHit *hit,
0223                                                  const GlobalPoint &gpos,
0224                                                  const GlobalVector &gdir,
0225                                                  bool mustProject) const {
0226   const TrackerSingleRecHit *stripHit = nullptr;
0227   if (typeid(*hit) == typeid(SiStripMatchedRecHit2D)) {
0228     const SiStripMatchedRecHit2D &mhit = static_cast<const SiStripMatchedRecHit2D &>(*hit);
0229     SiStripRecHit2D mono = mhit.monoHit();
0230     SiStripRecHit2D stereo = mhit.stereoHit();
0231     return testLastHit(&mono, gpos, gdir, true) && testLastHit(&stereo, gpos, gdir, true);
0232   } else if (typeid(*hit) == typeid(ProjectedSiStripRecHit2D)) {
0233     const ProjectedSiStripRecHit2D &mhit = static_cast<const ProjectedSiStripRecHit2D &>(*hit);
0234     const SiStripRecHit2D &orig = mhit.originalHit();
0235     return testLastHit(&orig, gpos, gdir, true);
0236   } else if ((stripHit = dynamic_cast<const TrackerSingleRecHit *>(hit)) != nullptr) {
0237     DetId detId = hit->geographicalId();
0238
0239     if (layerMask_[detId.subdetId()][0] == 0) {
0240       return true;  // no filtering here
0241     } else if (layerMask_[detId.subdetId()][0] == 2) {
0242       unsigned int ilayer = theTopology->layer(detId);
0243       if (layerMask_[detId.subdetId()][ilayer] == 0) {
0244         return true;  // no filtering here
0245       }
0246     }
0247
0248     const GeomDet *det = theTracker->idToDet(detId);
0249     LocalVector ldir = det->toLocal(gdir);
0250     LocalPoint lpos = det->toLocal(gpos);
0251     if (mustProject) {
0252       lpos -= ldir * lpos.z() / ldir.z();
0253     }
0254     int hitStrips;
0255     float hitPredPos;
0256     const SiStripCluster &cluster = stripHit->stripCluster();
0257     bool usable = theFilter->getSizes(detId, cluster, lpos, ldir, hitStrips, hitPredPos);
0258     if (!usable)
0259       return true;
0260
0261     INC_COUNTER(called_)
0262     const auto &ampls = cluster.amplitudes();
0263
0264     // pass-through of trivial case
0265     if (std::abs(hitPredPos) < 1.5f && hitStrips <= 2) {
0266       return true;
0267     }
0268
0269     if (!cluster.isFromApprox()) {
0270       // compute number of consecutive saturated strips
0271       // (i.e. with adc count >= 254, see SiStripCluster class for documentation)
0272       unsigned int thisSat = (ampls[0] >= 254), maxSat = thisSat;
0273       for (unsigned int i = 1, n = ampls.size(); i < n; ++i) {
0274         if (ampls[i] >= 254) {
0275           thisSat++;
0276         } else if (thisSat > 0) {
0277           maxSat = std::max<int>(maxSat, thisSat);
0278           thisSat = 0;
0279         }
0280       }
0281       if (thisSat > 0) {
0282         maxSat = std::max<int>(maxSat, thisSat);
0283       }
0284       if (maxSat >= maxNSat_) {
0285         INC_COUNTER(saturated_)
0286         return true;
0287       }
0288
0289       // trimming
0290       INC_COUNTER(test_)
0291       unsigned int hitStripsTrim = ampls.size();
0292       int sum = std::accumulate(ampls.begin(), ampls.end(), 0);
0293       uint8_t trimCut = std::min<uint8_t>(trimMaxADC_, std::floor(trimMaxFracTotal_ * sum));
0294       auto begin = ampls.begin();
0295       auto last = ampls.end() - 1;
0296       while (hitStripsTrim > 1 && (*begin < std::max<uint8_t>(trimCut, trimMaxFracNeigh_ * (*(begin + 1))))) {
0297         hitStripsTrim--;
0298         ++begin;
0299       }
0300       while (hitStripsTrim > 1 && (*last < std::max<uint8_t>(trimCut, trimMaxFracNeigh_ * (*(last - 1))))) {
0301         hitStripsTrim--;
0302         --last;
0303       }
0304
0305       if (hitStripsTrim < std::floor(std::abs(hitPredPos) - maxTrimmedSizeDiffNeg_)) {
0306         INC_COUNTER(failTooNarrow_)
0307         return false;
0308       } else if (hitStripsTrim <= std::ceil(std::abs(hitPredPos) + maxTrimmedSizeDiffPos_)) {
0309         INC_COUNTER(passTrim_)
0310         return true;
0311       }
0312
0313       const StripGeomDetUnit *stripDetUnit = dynamic_cast<const StripGeomDetUnit *>(det);
0314       if (det == nullptr) {
0315         edm::LogError("Strip not a StripGeomDetUnit?") << " on " << detId.rawId() << "\n";
0316         return true;
0317       }
0318
0319       float mip = 3.9 / (sistrip::MeVperADCStrip /
0320                          stripDetUnit->surface().bounds().thickness());  // 3.9 MeV/cm = ionization in silicon
0321       float mipnorm = mip / std::abs(ldir.z());
0322       ::SlidingPeakFinder pf(std::max<int>(2, std::ceil(std::abs(hitPredPos) + subclusterWindow_)));
0323       ::PeakFinderTest test(mipnorm,
0324                             detId(),
0325                             cluster.firstStrip(),
0326                             &*theNoise,
0327                             seedCutMIPs_,
0328                             seedCutSN_,
0329                             subclusterCutMIPs_,
0330                             subclusterCutSN_);
0331       if (pf.apply(cluster.amplitudes(), test)) {
0332         INC_COUNTER(passSC_)
0333         return true;
0334       } else {
0335         INC_COUNTER(failTooLarge_)
0336         return false;
0337       }
0338     } else {
0339       return cluster.filter();
0340     }
0341   }
0342   return true;
0343 }
0344
0345 void StripSubClusterShapeFilterBase::setEventBase(const edm::Event &event, const edm::EventSetup &es) {
0346   // Get tracker geometry
0347   theTracker = es.getHandle(geomToken_);
0348   theFilter = es.getHandle(csfToken_);
0349
0350   //Retrieve tracker topology from geometry
0351   theTopology = es.getHandle(topoToken_);
0352   theNoise = es.getHandle(stripNoiseToken_);
0353 }
0354
0355 /*****************************************************************************/
0356 bool StripSubClusterShapeTrajectoryFilter::toBeContinued(Trajectory &trajectory) const {
0357   throw cms::Exception("toBeContinued(Traj) instead of toBeContinued(TempTraj)");
0358 }
0359
0360 /*****************************************************************************/
0361 bool StripSubClusterShapeTrajectoryFilter::testLastHit(const TrajectoryMeasurement &last) const {
0362   const TrackingRecHit *hit = last.recHit()->hit();
0363   if (!last.updatedState().isValid())
0364     return true;
0365   if (hit == nullptr || !hit->isValid())
0366     return true;
0367   if (hit->geographicalId().subdetId() < SiStripDetId::TIB)
0368     return true;  // we look only at strips for now
0369   return testLastHit(hit, last.updatedState(), false);
0370 }
0371
0372 /*****************************************************************************/
0373 bool StripSubClusterShapeTrajectoryFilter::toBeContinued(TempTrajectory &trajectory) const {
0374   const TempTrajectory::DataContainer &tms = trajectory.measurements();
0375   return testLastHit(*tms.rbegin());
0376 }
0377
0378 /*****************************************************************************/
0379 bool StripSubClusterShapeTrajectoryFilter::qualityFilter(const Trajectory &trajectory) const {
0380   const Trajectory::DataContainer &tms = trajectory.measurements();
0381   return testLastHit(*tms.rbegin());
0382 }
0383
0384 /*****************************************************************************/
0385 bool StripSubClusterShapeTrajectoryFilter::qualityFilter(const TempTrajectory &trajectory) const {
0386   const TempTrajectory::DataContainer &tms = trajectory.measurements();
0387   return testLastHit(*tms.rbegin());
0388 }
0389
0390 /*****************************************************************************/
0391 StripSubClusterShapeSeedFilter::StripSubClusterShapeSeedFilter(const edm::ParameterSet &iConfig,
0392                                                                edm::ConsumesCollector &iC)
0393     : StripSubClusterShapeFilterBase(iConfig, iC),
0394       filterAtHelixStage_(iConfig.getParameter<bool>("FilterAtHelixStage"))
0395
0396 {
0397   if (filterAtHelixStage_)
0398     edm::LogError("Configuration")
0399         << "StripSubClusterShapeSeedFilter: FilterAtHelixStage is not yet working correctly.\n";
0400 }
0401
0402 /*****************************************************************************/
0403 bool StripSubClusterShapeSeedFilter::compatible(const TrajectoryStateOnSurface &tsos,
0404                                                 SeedingHitSet::ConstRecHitPointer thit) const {
0405   if (filterAtHelixStage_)
0406     return true;
0407   const TrackingRecHit *hit = thit->hit();
0408   if (hit == nullptr || !hit->isValid())
0409     return true;
0410   if (hit->geographicalId().subdetId() < SiStripDetId::TIB)
0411     return true;  // we look only at strips for now
0412   return testLastHit(hit, tsos, false);
0413 }
0414
0415 bool StripSubClusterShapeSeedFilter::compatible(const SeedingHitSet &hits,
0416                                                 const GlobalTrajectoryParameters &helixStateAtVertex,
0417                                                 const FastHelix &helix) const {
0418   if (!filterAtHelixStage_)
0419     return true;
0420
0421   if (!helix.isValid()              //check still if it's a straight line, which are OK
0422       && !helix.circle().isLine())  //complain if it's not even a straight line
0423     edm::LogWarning("InvalidHelix") << "StripSubClusterShapeSeedFilter helix is not valid, result is bad";
0424
0425   float xc = helix.circle().x0(), yc = helix.circle().y0();
0426
0427   GlobalPoint vertex = helixStateAtVertex.position();
0428   GlobalVector momvtx = helixStateAtVertex.momentum();
0429   float x0 = vertex.x(), y0 = vertex.y();
0430   for (unsigned int i = 0, n = hits.size(); i < n; ++i) {
0431     auto const &hit = *hits[i];
0432     if (hit.geographicalId().subdetId() < SiStripDetId::TIB)
0433       continue;
0434
0435     GlobalPoint pos = hit.globalPosition();
0436     float x1 = pos.x(), y1 = pos.y(), dx1 = x1 - xc, dy1 = y1 - yc;
0437
0438     // now figure out the proper tangent vector
0439     float perpx = -dy1, perpy = dx1;
0440     if (perpx * (x1 - x0) + perpy * (y1 - y0) < 0) {
0441       perpy = -perpy;
0442       perpx = -perpx;
0443     }
0444
0445     // now normalize (perpx, perpy, 1.0) to unity
0446     float pnorm = 1.0 / std::sqrt(perpx * perpx + perpy * perpy + 1);
0447     GlobalVector gdir(perpx * pnorm, perpy * pnorm, (momvtx.z() > 0 ? pnorm : -pnorm));
0448
0449     if (!testLastHit(&hit, pos, gdir)) {
0450       return false;  // not yet
0451     }
0452   }
0453   return true;
0454 }