File indexing completed on 2023-03-17 11:22:46
0001 #include "RecoTracker/SpecialSeedGenerators/interface/GenericTripletGenerator.h"
0002 #include "RecoTracker/TkSeedGenerator/interface/FastCircle.h"
0003 typedef SeedingHitSet::ConstRecHitPointer SeedingHit;
0004
0005 #include "FWCore/Framework/interface/ConsumesCollector.h"
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 #include <map>
0008
0009 GenericTripletGenerator::GenericTripletGenerator(const edm::ParameterSet& conf, edm::ConsumesCollector& iC)
0010 : theSeedingLayerToken(iC.consumes<SeedingLayerSetsHits>(conf.getParameter<edm::InputTag>("LayerSrc"))) {
0011 edm::LogInfo("CtfSpecialSeedGenerator|GenericTripletGenerator") << "Constructing GenericTripletGenerator";
0012 }
0013
0014 const OrderedSeedingHits& GenericTripletGenerator::run(const TrackingRegion& region,
0015 const edm::Event& e,
0016 const edm::EventSetup& es) {
0017 hitTriplets.clear();
0018 hitTriplets.reserve(0);
0019 edm::Handle<SeedingLayerSetsHits> hlayers;
0020 e.getByToken(theSeedingLayerToken, hlayers);
0021 const SeedingLayerSetsHits& layers = *hlayers;
0022 if (layers.numberOfLayersInSet() != 3)
0023 throw cms::Exception("CtfSpecialSeedGenerator")
0024 << "You are using " << layers.numberOfLayersInSet() << " layers in set instead of 3 ";
0025 std::map<float, OrderedHitTriplet> radius_triplet_map;
0026 for (SeedingLayerSetsHits::SeedingLayerSet ls : layers) {
0027 auto innerHits = region.hits(ls[0]);
0028
0029 auto middleHits = region.hits(ls[1]);
0030
0031 auto outerHits = region.hits(ls[2]);
0032
0033
0034 for (auto iOuterHit = outerHits.begin(); iOuterHit != outerHits.end(); iOuterHit++) {
0035 for (auto iMiddleHit = middleHits.begin(); iMiddleHit != middleHits.end(); iMiddleHit++) {
0036 for (auto iInnerHit = innerHits.begin(); iInnerHit != innerHits.end(); iInnerHit++) {
0037
0038
0039
0040
0041
0042
0043
0044 OrderedHitTriplet oht(&(**iInnerHit), &(**iMiddleHit), &(**iOuterHit));
0045 std::pair<bool, float> val_radius = qualityFilter(oht, radius_triplet_map, ls);
0046 if (val_radius.first) {
0047
0048
0049
0050
0051 hitTriplets.push_back(oht);
0052 radius_triplet_map.insert(std::make_pair(val_radius.second, oht));
0053 }
0054 }
0055 }
0056 }
0057 }
0058
0059 return hitTriplets;
0060 }
0061
0062 std::pair<bool, float> GenericTripletGenerator::qualityFilter(const OrderedHitTriplet& oht,
0063 const std::map<float, OrderedHitTriplet>& map,
0064 const SeedingLayerSetsHits::SeedingLayerSet& ls) const {
0065
0066 GlobalPoint innerpos = oht.inner()->globalPosition();
0067 GlobalPoint middlepos = oht.middle()->globalPosition();
0068 GlobalPoint outerpos = oht.outer()->globalPosition();
0069 std::vector<const TrackingRecHit*> ohttrh;
0070 ohttrh.push_back(&(*(oht.inner())));
0071 ohttrh.push_back(&(*(oht.middle())));
0072 ohttrh.push_back(&(*(oht.outer())));
0073 std::vector<const TrackingRecHit*>::const_iterator ioht;
0074
0075
0076 float limit_phi_distance1 = sqrt((middlepos.x() - outerpos.x()) * (middlepos.x() - outerpos.x()) +
0077 (middlepos.y() - outerpos.y()) * (middlepos.y() - outerpos.y())) /
0078 middlepos.mag();
0079 float limit_phi_distance2 = sqrt((middlepos.x() - innerpos.x()) * (middlepos.x() - innerpos.x()) +
0080 (middlepos.y() - innerpos.y()) * (middlepos.y() - innerpos.y())) /
0081 innerpos.mag();
0082
0083
0084 if (fabs(outerpos.phi() - middlepos.phi()) > fabs(atan(limit_phi_distance1)) ||
0085 fabs(innerpos.phi() - middlepos.phi()) > fabs(atan(limit_phi_distance2))) {
0086
0087 return std::make_pair(false, 0.);
0088 }
0089
0090
0091
0092
0093
0094
0095
0096
0097
0098
0099
0100
0101
0102
0103
0104
0105
0106
0107 FastCircle circle(innerpos, middlepos, outerpos);
0108 if (circle.rho() < 200 && circle.rho() != 0)
0109 return std::make_pair(false, circle.rho());
0110
0111
0112 std::map<float, OrderedHitTriplet>::const_iterator lower_bound = map.lower_bound((1 - 0.01) * circle.rho());
0113 std::map<float, OrderedHitTriplet>::const_iterator upper_bound = map.upper_bound((1 + 0.01) * circle.rho());
0114 std::map<float, OrderedHitTriplet>::const_iterator iter;
0115 for (iter = lower_bound; iter != upper_bound && iter->first <= upper_bound->first; iter++) {
0116 int shared = 0;
0117 std::vector<const TrackingRecHit*> curtrh;
0118 curtrh.push_back(&*(iter->second.inner()));
0119 curtrh.push_back(&*(iter->second.middle()));
0120 curtrh.push_back(&*(iter->second.outer()));
0121 std::vector<const TrackingRecHit*>::const_iterator curiter;
0122 for (curiter = curtrh.begin(); curiter != curtrh.end(); curiter++) {
0123 for (ioht = ohttrh.begin(); ioht != ohttrh.end(); ioht++) {
0124 if ((*ioht)->geographicalId() == (*curiter)->geographicalId() &&
0125 ((*ioht)->localPosition() - (*curiter)->localPosition()).mag() < 1e-5)
0126 shared++;
0127 }
0128 }
0129 if (shared > 1)
0130 return std::make_pair(false, circle.rho());
0131 }
0132
0133 return std::make_pair(true, circle.rho());
0134 }