Back to home page

Project CMSSW displayed by LXR
 
 

    

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 &ampls, 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 }  // namespace
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;  // no filtering here
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;  // no filtering here
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 &ampls = cluster.amplitudes();
0262 
0263     // pass-through of trivial case
0264     if (std::abs(hitPredPos) < 1.5f && hitStrips <= 2) {
0265       return true;
0266     }
0267 
0268     // compute number of consecutive saturated strips
0269     // (i.e. with adc count >= 254, see SiStripCluster class for documentation)
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     // trimming
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;  // conversion constant from ADC counts to MeV for the SiStrip detector
0318     float mip =
0319         3.9 / (MeVperADCStrip / stripDetUnit->surface().bounds().thickness());  // 3.9 MeV/cm = ionization in silicon
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   // Get tracker geometry
0343   theTracker = es.getHandle(geomToken_);
0344   theFilter = es.getHandle(csfToken_);
0345 
0346   //Retrieve tracker topology from geometry
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;  // we look only at strips for now
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;  // we look only at strips for now
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()              //check still if it's a straight line, which are OK
0418       && !helix.circle().isLine())  //complain if it's not even a straight line
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     // now figure out the proper tangent vector
0435     float perpx = -dy1, perpy = dx1;
0436     if (perpx * (x1 - x0) + perpy * (y1 - y0) < 0) {
0437       perpy = -perpy;
0438       perpx = -perpx;
0439     }
0440 
0441     // now normalize (perpx, perpy, 1.0) to unity
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;  // not yet
0447     }
0448   }
0449   return true;
0450 }
[ 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!