File indexing completed on 2024-04-06 12:28:55
0001 #include "HitExtractorSTRP.h"
0002 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0003
0004 #include "DataFormats/Common/interface/Handle.h"
0005 #include "FWCore/Framework/interface/Event.h"
0006
0007 #include "FWCore/Framework/interface/EventSetup.h"
0008
0009 #include "DataFormats/Common/interface/ContainerMask.h"
0010
0011 #include "TrackingTools/TransientTrackingRecHit/interface/TrackingRecHitProjector.h"
0012 #include "RecoTracker/TransientTrackingRecHit/interface/ProjectedRecHit2D.h"
0013 #include "DataFormats/SiStripCluster/interface/SiStripClusterTools.h"
0014
0015 #include <tuple>
0016
0017 #include <iostream>
0018
0019 using namespace ctfseeding;
0020 using namespace std;
0021 using namespace edm;
0022
0023 HitExtractorSTRP::HitExtractorSTRP(GeomDetEnumerators::SubDetector subdet,
0024 TrackerDetSide side,
0025 int idLayer,
0026 float iminGoodCharge,
0027 edm::ConsumesCollector& iC)
0028 : theLayerSubDet(subdet),
0029 theSide(side),
0030 theIdLayer(idLayer),
0031 minAbsZ(0),
0032 theMinRing(1),
0033 theMaxRing(0),
0034 theTtopo(iC.esConsumes()),
0035 hasMatchedHits(false),
0036 hasRPhiHits(false),
0037 hasStereoHits(false),
0038 hasVectorHits(false),
0039 hasRingSelector(false),
0040 hasSimpleRphiHitsCleaner(true) {
0041 minGoodCharge = iminGoodCharge;
0042 if (minGoodCharge > 0)
0043 skipClusters = true;
0044 }
0045
0046 void HitExtractorSTRP::useSkipClusters_(const edm::InputTag& m, edm::ConsumesCollector& iC) {
0047 theSkipClusters = iC.consumes<SkipClustersCollection>(m);
0048 theSkipPhase2Clusters = iC.consumes<SkipPhase2ClustersCollection>(m);
0049 }
0050
0051 void HitExtractorSTRP::useRingSelector(int minRing, int maxRing) {
0052 hasRingSelector = true;
0053 theMinRing = minRing;
0054 theMaxRing = maxRing;
0055 }
0056
0057 bool HitExtractorSTRP::ringRange(int ring) const {
0058 if (!hasRingSelector)
0059 return true;
0060 return (ring >= theMinRing) & (ring <= theMaxRing);
0061 }
0062
0063 bool HitExtractorSTRP::skipThis(
0064 DetId id,
0065 OmniClusterRef const& clus,
0066 edm::Handle<edm::ContainerMask<edmNew::DetSetVector<SiStripCluster> > >& stripClusterMask) const {
0067 if (maskCluster && (stripClusterMask->mask(clus.key())))
0068 return true;
0069
0070 if UNLIKELY (minGoodCharge <= 0)
0071 return false;
0072 return siStripClusterTools::chargePerCM(id, *clus.cluster_strip()) <= minGoodCharge;
0073 }
0074
0075 std::pair<bool, ProjectedSiStripRecHit2D*> HitExtractorSTRP::skipThis(
0076 const TkTransientTrackingRecHitBuilder& ttrhBuilder,
0077 TkHitRef matched,
0078 edm::Handle<edm::ContainerMask<edmNew::DetSetVector<SiStripCluster> > >& stripClusterMask) const {
0079 const SiStripMatchedRecHit2D& hit = (SiStripMatchedRecHit2D const&)(matched);
0080
0081 assert(dynamic_cast<SiStripMatchedRecHit2D const*>(&matched));
0082
0083 auto id = hit.geographicalId();
0084 ProjectedSiStripRecHit2D* replaceMe = nullptr;
0085 bool rejectSt = skipThis(id, hit.stereoClusterRef(), stripClusterMask);
0086 bool rejectMono = skipThis(id, hit.monoClusterRef(), stripClusterMask);
0087
0088 if ((!rejectSt) & (!rejectMono)) {
0089
0090 return std::make_pair(false, replaceMe);
0091 }
0092
0093 if (failProjection || (rejectSt & rejectMono)) {
0094
0095 return std::make_pair(true, replaceMe);
0096 }
0097
0098
0099
0100 auto cloner = ttrhBuilder.cloner();
0101 replaceMe = cloner.project(hit, rejectSt, TrajectoryStateOnSurface()).release();
0102 if (rejectSt)
0103 LogDebug("HitExtractorSTRP") << "a matched hit is partially masked, and the mono hit got projected onto: "
0104 << replaceMe->geographicalId().rawId() << " key: " << hit.monoClusterRef().key();
0105 else
0106 LogDebug("HitExtractorSTRP") << "a matched hit is partially masked, and the stereo hit got projected onto: "
0107 << replaceMe->geographicalId().rawId() << " key: " << hit.stereoClusterRef().key();
0108
0109 return std::make_pair(true, replaceMe);
0110 }
0111
0112 void HitExtractorSTRP::cleanedOfClusters(const TkTransientTrackingRecHitBuilder& ttrhBuilder,
0113 const edm::Event& ev,
0114 HitExtractor::Hits& hits,
0115 bool matched,
0116 unsigned int cleanFrom) const {
0117 unsigned int skipped = 0;
0118 unsigned int projected = 0;
0119 if (hasMatchedHits || hasRPhiHits || hasStereoHits) {
0120 LogTrace("HitExtractorSTRP") << "getting " << hits.size() << " strip hit in input.";
0121 edm::Handle<SkipClustersCollection> stripClusterMask;
0122 if (maskCluster)
0123 ev.getByToken(theSkipClusters, stripClusterMask);
0124 for (unsigned int iH = cleanFrom; iH < hits.size(); ++iH) {
0125 assert(hits[iH]->isValid());
0126 auto id = hits[iH]->geographicalId();
0127 if (matched) {
0128 auto [replace, replaceMe] = skipThis(ttrhBuilder, *hits[iH], stripClusterMask);
0129 if (replace) {
0130 if (!replaceMe) {
0131 LogTrace("HitExtractorSTRP") << "skipping a matched hit on :" << hits[iH]->geographicalId().rawId();
0132 skipped++;
0133 } else {
0134 projected++;
0135 }
0136 hits[iH].reset(replaceMe);
0137 if (replaceMe == nullptr)
0138 assert(hits[iH].empty());
0139 else
0140 assert(hits[iH].isOwn());
0141 }
0142 } else if (skipThis(id, hits[iH]->firstClusterRef(), stripClusterMask)) {
0143 LogTrace("HitExtractorSTRP") << "skipping a hit on :" << hits[iH]->geographicalId().rawId() << " key: ";
0144 skipped++;
0145 hits[iH].reset();
0146 }
0147 }
0148 }
0149 if (hasVectorHits) {
0150 LogTrace("HitExtractorSTRP") << "getting " << hits.size() << " vector hit in input.";
0151 edm::Handle<SkipPhase2ClustersCollection> ph2ClusterMask;
0152 if (maskCluster)
0153 ev.getByToken(theSkipPhase2Clusters, ph2ClusterMask);
0154 for (unsigned int iH = cleanFrom; iH < hits.size(); ++iH) {
0155 LogTrace("HitExtractorSTRP") << "analizing hit on :" << hits[iH]->geographicalId().rawId();
0156 assert(hits[iH]->isValid());
0157 const VectorHit& vhit = dynamic_cast<VectorHit const&>(*hits[iH]);
0158 LogTrace("HitExtractorSTRP") << " key lower: " << vhit.lowerClusterRef().key()
0159 << " and key upper: " << vhit.upperClusterRef().key();
0160 LogTrace("HitExtractorSTRP") << " key lower: " << hits[iH]->firstClusterRef().key();
0161
0162
0163 if (maskCluster &&
0164 (ph2ClusterMask->mask(vhit.lowerClusterRef().key()) || ph2ClusterMask->mask(vhit.upperClusterRef().key()))) {
0165 LogTrace("HitExtractorSTRP") << "skipping a vector hit on :" << hits[iH]->geographicalId().rawId()
0166 << " key lower: " << vhit.lowerClusterRef().key()
0167 << " and key upper: " << vhit.upperClusterRef().key();
0168 skipped++;
0169 hits[iH].reset();
0170 }
0171 }
0172 }
0173
0174
0175 auto last = std::remove_if(hits.begin() + cleanFrom, hits.end(), [](HitPointer const& p) { return p.empty(); });
0176 hits.resize(last - hits.begin());
0177
0178 LogTrace("HitExtractorSTRP") << "skipped :" << skipped << " rechits because of clusters and projected: " << projected;
0179 }
0180
0181 HitExtractor::Hits HitExtractorSTRP::hits(const TkTransientTrackingRecHitBuilder& ttrhBuilder,
0182 const edm::Event& ev,
0183 const edm::EventSetup& es) const {
0184 LogDebug("HitExtractorSTRP") << "HitExtractorSTRP::hits";
0185 HitExtractor::Hits result;
0186 unsigned int cleanFrom = 0;
0187
0188
0189 const TrackerTopology* const tTopo = &es.getData(theTtopo);
0190
0191
0192
0193
0194 if (theLayerSubDet == GeomDetEnumerators::TIB) {
0195 LogTrace("HitExtractorSTRP") << "Getting hits into the TIB";
0196 if (hasMatchedHits) {
0197 edm::Handle<SiStripMatchedRecHit2DCollection> matchedHits;
0198 ev.getByToken(theMatchedHits, matchedHits);
0199 if (skipClusters)
0200 cleanFrom = result.size();
0201 range2SeedingHits(*matchedHits, result, tTopo->tibDetIdLayerComparator(theIdLayer));
0202 if (skipClusters)
0203 cleanedOfClusters(ttrhBuilder, ev, result, true, cleanFrom);
0204 }
0205 if (hasRPhiHits) {
0206 edm::Handle<SiStripRecHit2DCollection> rphiHits;
0207 ev.getByToken(theRPhiHits, rphiHits);
0208 if (hasMatchedHits) {
0209 if (!hasSimpleRphiHitsCleaner) {
0210 if (skipClusters)
0211 cleanFrom = result.size();
0212 range2SeedingHits(*rphiHits, result, tTopo->tibDetIdLayerComparator(theIdLayer));
0213 if (skipClusters)
0214 cleanedOfClusters(ttrhBuilder, ev, result, false, cleanFrom);
0215 }
0216 } else {
0217 if (skipClusters)
0218 cleanFrom = result.size();
0219 range2SeedingHits(*rphiHits, result, tTopo->tibDetIdLayerComparator(theIdLayer));
0220 if (skipClusters)
0221 cleanedOfClusters(ttrhBuilder, ev, result, false, cleanFrom);
0222 }
0223 }
0224 if (hasStereoHits) {
0225 edm::Handle<SiStripRecHit2DCollection> stereoHits;
0226 ev.getByToken(theStereoHits, stereoHits);
0227 if (skipClusters)
0228 cleanFrom = result.size();
0229 range2SeedingHits(*stereoHits, result, tTopo->tibDetIdLayerComparator(theIdLayer));
0230 if (skipClusters)
0231 cleanedOfClusters(ttrhBuilder, ev, result, false, cleanFrom);
0232 }
0233 if (hasVectorHits) {
0234 LogError("HitExtractorSTRP") << "TIB is not supposed to be in Phase2 TRK detector configuration. What follows "
0235 "have never been checked before! ";
0236 auto const& vectorHits = ev.get(theVectorHits);
0237 if (skipClusters)
0238 cleanFrom = result.size();
0239 range2SeedingHits(vectorHits, result, tTopo->tibDetIdLayerComparator(theIdLayer));
0240 if (skipClusters)
0241 cleanedOfClusters(ttrhBuilder, ev, result, false, cleanFrom);
0242 }
0243
0244 }
0245
0246
0247
0248
0249 else if (theLayerSubDet == GeomDetEnumerators::TID) {
0250 LogTrace("HitExtractorSTRP") << "Getting hits into the TID";
0251 if (hasMatchedHits) {
0252 edm::Handle<SiStripMatchedRecHit2DCollection> matchedHits;
0253 ev.getByToken(theMatchedHits, matchedHits);
0254 if (skipClusters)
0255 cleanFrom = result.size();
0256 auto getter = tTopo->tidDetIdWheelComparator(static_cast<unsigned int>(theSide), theIdLayer);
0257 SiStripMatchedRecHit2DCollection::Range range = matchedHits->equal_range(getter.first, getter.second);
0258 for (SiStripMatchedRecHit2DCollection::const_iterator it = range.first; it != range.second; ++it) {
0259 int ring = tTopo->tidRing(it->detId());
0260 if (!ringRange(ring))
0261 continue;
0262 for (SiStripMatchedRecHit2DCollection::DetSet::const_iterator hit = it->begin(), end = it->end(); hit != end;
0263 ++hit) {
0264 result.emplace_back(*hit);
0265 }
0266 }
0267 if (skipClusters)
0268 cleanedOfClusters(ttrhBuilder, ev, result, true, cleanFrom);
0269 }
0270 if (hasRPhiHits) {
0271 edm::Handle<SiStripRecHit2DCollection> rphiHits;
0272 ev.getByToken(theRPhiHits, rphiHits);
0273 if (skipClusters)
0274 cleanFrom = result.size();
0275 auto getter = tTopo->tidDetIdWheelComparator(static_cast<unsigned int>(theSide), theIdLayer);
0276 SiStripRecHit2DCollection::Range range = rphiHits->equal_range(getter.first, getter.second);
0277 for (SiStripRecHit2DCollection::const_iterator it = range.first; it != range.second; ++it) {
0278 int ring = tTopo->tidRing(it->detId());
0279 if (!ringRange(ring))
0280 continue;
0281 if ((SiStripDetId(it->detId()).partnerDetId() != 0) && hasSimpleRphiHitsCleaner)
0282 continue;
0283 for (SiStripRecHit2DCollection::DetSet::const_iterator hit = it->begin(), end = it->end(); hit != end; ++hit) {
0284 result.emplace_back(*hit);
0285 }
0286 }
0287 if (skipClusters)
0288 cleanedOfClusters(ttrhBuilder, ev, result, false, cleanFrom);
0289 }
0290 if (hasStereoHits) {
0291 edm::Handle<SiStripRecHit2DCollection> stereoHits;
0292 ev.getByToken(theStereoHits, stereoHits);
0293 if (skipClusters)
0294 cleanFrom = result.size();
0295 auto getter = tTopo->tidDetIdWheelComparator(static_cast<unsigned int>(theSide), theIdLayer);
0296 SiStripRecHit2DCollection::Range range = stereoHits->equal_range(getter.first, getter.second);
0297 for (SiStripRecHit2DCollection::const_iterator it = range.first; it != range.second; ++it) {
0298 int ring = tTopo->tidRing(it->detId());
0299 if (!ringRange(ring))
0300 continue;
0301 for (SiStripRecHit2DCollection::DetSet::const_iterator hit = it->begin(), end = it->end(); hit != end; ++hit) {
0302 result.emplace_back(*hit);
0303 }
0304 }
0305 if (skipClusters)
0306 cleanedOfClusters(ttrhBuilder, ev, result, false, cleanFrom);
0307 }
0308 if (hasVectorHits) {
0309 LogTrace("HitExtractorSTRP") << "Getting vector hits for IdLayer " << theIdLayer;
0310 auto const& vectorHits = ev.get(theVectorHits);
0311
0312 if (skipClusters)
0313 cleanFrom = result.size();
0314 auto getter = tTopo->tidDetIdWheelComparator(static_cast<unsigned int>(theSide), theIdLayer);
0315 VectorHitCollection::Range range = vectorHits.equal_range(getter.first, getter.second);
0316 for (VectorHitCollection::const_iterator it = range.first; it != range.second; ++it) {
0317 int ring = tTopo->tidRing(it->detId());
0318 if (!ringRange(ring))
0319 continue;
0320 for (VectorHitCollection::DetSet::const_iterator hit = it->begin(), end = it->end(); hit != end; ++hit) {
0321 result.emplace_back(*hit);
0322 }
0323 }
0324 LogTrace("HitExtractorSTRP") << "result size value:" << result.size();
0325 if (skipClusters)
0326 cleanedOfClusters(ttrhBuilder, ev, result, false, cleanFrom);
0327 }
0328 }
0329
0330
0331
0332 else if (theLayerSubDet == GeomDetEnumerators::TOB) {
0333 LogTrace("HitExtractorSTRP") << "Getting hits into the TOB";
0334 if (hasMatchedHits) {
0335 edm::Handle<SiStripMatchedRecHit2DCollection> matchedHits;
0336 ev.getByToken(theMatchedHits, matchedHits);
0337 if (skipClusters)
0338 cleanFrom = result.size();
0339 if (minAbsZ > 0.) {
0340 auto getter = tTopo->tobDetIdLayerComparator(theIdLayer);
0341 SiStripMatchedRecHit2DCollection::Range range = matchedHits->equal_range(getter.first, getter.second);
0342 for (SiStripMatchedRecHit2DCollection::const_iterator it = range.first; it != range.second; ++it) {
0343 for (SiStripMatchedRecHit2DCollection::DetSet::const_iterator hit = it->begin(), end = it->end(); hit != end;
0344 ++hit) {
0345 if (fabs(hit->globalPosition().z()) >= minAbsZ)
0346 result.emplace_back(*hit);
0347 }
0348 }
0349 } else {
0350 range2SeedingHits(*matchedHits, result, tTopo->tobDetIdLayerComparator(theIdLayer));
0351 }
0352 if (skipClusters)
0353 cleanedOfClusters(ttrhBuilder, ev, result, true, cleanFrom);
0354 }
0355 if (hasRPhiHits) {
0356 edm::Handle<SiStripRecHit2DCollection> rphiHits;
0357 ev.getByToken(theRPhiHits, rphiHits);
0358 if (hasMatchedHits) {
0359 if (!hasSimpleRphiHitsCleaner) {
0360 if (skipClusters)
0361 cleanFrom = result.size();
0362 range2SeedingHits(*rphiHits, result, tTopo->tobDetIdLayerComparator(theIdLayer));
0363 if (skipClusters)
0364 cleanedOfClusters(ttrhBuilder, ev, result, false, cleanFrom);
0365 }
0366 } else {
0367 if (skipClusters)
0368 cleanFrom = result.size();
0369 range2SeedingHits(*rphiHits, result, tTopo->tobDetIdLayerComparator(theIdLayer));
0370 if (skipClusters)
0371 cleanedOfClusters(ttrhBuilder, ev, result, false, cleanFrom);
0372 }
0373 }
0374 if (hasStereoHits) {
0375 edm::Handle<SiStripRecHit2DCollection> stereoHits;
0376 ev.getByToken(theStereoHits, stereoHits);
0377 if (skipClusters)
0378 cleanFrom = result.size();
0379 range2SeedingHits(*stereoHits, result, tTopo->tobDetIdLayerComparator(theIdLayer));
0380 if (skipClusters)
0381 cleanedOfClusters(ttrhBuilder, ev, result, false, cleanFrom);
0382 }
0383 if (hasVectorHits) {
0384 LogTrace("HitExtractorSTRP") << "Getting vector hits for IdLayer " << theIdLayer;
0385 edm::Handle<VectorHitCollection> vectorHits;
0386 ev.getByToken(theVectorHits, vectorHits);
0387
0388 if (skipClusters)
0389 cleanFrom = result.size();
0390 range2SeedingHits(*vectorHits, result, tTopo->tobDetIdLayerComparator(theIdLayer));
0391 if (skipClusters)
0392 cleanedOfClusters(ttrhBuilder, ev, result, false, cleanFrom);
0393 }
0394
0395 }
0396
0397
0398
0399
0400 else if (theLayerSubDet == GeomDetEnumerators::TEC) {
0401 LogTrace("HitExtractorSTRP") << "Getting hits into the TEC";
0402 if (hasMatchedHits) {
0403 edm::Handle<SiStripMatchedRecHit2DCollection> matchedHits;
0404 ev.getByToken(theMatchedHits, matchedHits);
0405 if (skipClusters)
0406 cleanFrom = result.size();
0407 auto getter = tTopo->tecDetIdWheelComparator(static_cast<unsigned int>(theSide), theIdLayer);
0408 SiStripMatchedRecHit2DCollection::Range range = matchedHits->equal_range(getter.first, getter.second);
0409 for (SiStripMatchedRecHit2DCollection::const_iterator it = range.first; it != range.second; ++it) {
0410 int ring = tTopo->tecRing(it->detId());
0411 if (!ringRange(ring))
0412 continue;
0413 for (SiStripMatchedRecHit2DCollection::DetSet::const_iterator hit = it->begin(), end = it->end(); hit != end;
0414 ++hit) {
0415 result.emplace_back(*hit);
0416 }
0417 }
0418 if (skipClusters)
0419 cleanedOfClusters(ttrhBuilder, ev, result, true, cleanFrom);
0420 }
0421 if (hasRPhiHits) {
0422 edm::Handle<SiStripRecHit2DCollection> rphiHits;
0423 ev.getByToken(theRPhiHits, rphiHits);
0424 if (skipClusters)
0425 cleanFrom = result.size();
0426 auto getter = tTopo->tecDetIdWheelComparator(static_cast<unsigned int>(theSide), theIdLayer);
0427 SiStripRecHit2DCollection::Range range = rphiHits->equal_range(getter.first, getter.second);
0428 for (SiStripRecHit2DCollection::const_iterator it = range.first; it != range.second; ++it) {
0429 int ring = tTopo->tecRing(it->detId());
0430 if (!ringRange(ring))
0431 continue;
0432 if ((SiStripDetId(it->detId()).partnerDetId() != 0) && hasSimpleRphiHitsCleaner)
0433 continue;
0434 for (SiStripRecHit2DCollection::DetSet::const_iterator hit = it->begin(), end = it->end(); hit != end; ++hit) {
0435 result.emplace_back(*hit);
0436 }
0437 }
0438 if (skipClusters)
0439 cleanedOfClusters(ttrhBuilder, ev, result, false, cleanFrom);
0440 }
0441 if (hasStereoHits) {
0442 edm::Handle<SiStripRecHit2DCollection> stereoHits;
0443 ev.getByToken(theStereoHits, stereoHits);
0444 if (skipClusters)
0445 cleanFrom = result.size();
0446 auto getter = tTopo->tecDetIdWheelComparator(static_cast<unsigned int>(theSide), theIdLayer);
0447 SiStripRecHit2DCollection::Range range = stereoHits->equal_range(getter.first, getter.second);
0448 for (SiStripRecHit2DCollection::const_iterator it = range.first; it != range.second; ++it) {
0449 int ring = tTopo->tecRing(it->detId());
0450 if (!ringRange(ring))
0451 continue;
0452 for (SiStripRecHit2DCollection::DetSet::const_iterator hit = it->begin(), end = it->end(); hit != end; ++hit) {
0453 result.emplace_back(*hit);
0454 }
0455 }
0456 if (skipClusters)
0457 cleanedOfClusters(ttrhBuilder, ev, result, false, cleanFrom);
0458 }
0459 if (hasVectorHits) {
0460 LogError("HitExtractorSTRP") << "TEC is not supposed to be in Phase2 TRK detector configuration. What follows "
0461 "have never been checked before! ";
0462 edm::Handle<VectorHitCollection> vectorHits;
0463 ev.getByToken(theVectorHits, vectorHits);
0464 if (skipClusters)
0465 cleanFrom = result.size();
0466 auto getter = tTopo->tidDetIdWheelComparator(static_cast<unsigned int>(theSide), theIdLayer);
0467 VectorHitCollection::Range range = vectorHits->equal_range(getter.first, getter.second);
0468 for (VectorHitCollection::const_iterator it = range.first; it != range.second; ++it) {
0469 int ring = tTopo->tidRing(it->detId());
0470 if (!ringRange(ring))
0471 continue;
0472 for (VectorHitCollection::DetSet::const_iterator hit = it->begin(), end = it->end(); hit != end; ++hit) {
0473 result.emplace_back(*hit);
0474 }
0475 }
0476 if (skipClusters)
0477 cleanedOfClusters(ttrhBuilder, ev, result, false, cleanFrom);
0478 }
0479 }
0480
0481 LogDebug("HitExtractorSTRP") << " giving: " << result.size() << " out for charge cut " << minGoodCharge;
0482 return result;
0483 }