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