File indexing completed on 2023-03-17 11:23:02
0001 #include "RecoTracker/TkSeedingLayers/interface/SeedingLayerSetsBuilder.h"
0002 #include "FWCore/Framework/interface/Event.h"
0003 #include "FWCore/Framework/interface/EventSetup.h"
0004 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0005 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0006 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0007 #include "RecoTracker/TkDetLayers/interface/GeometricSearchTracker.h"
0008 #include "RecoTracker/Record/interface/TrackerRecoGeometryRecord.h"
0009
0010 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
0011 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h"
0012 #include "RecoTracker/TransientTrackingRecHit/interface/TkTransientTrackingRecHitBuilder.h"
0013
0014 #include "TrackingTools/DetLayers/interface/BarrelDetLayer.h"
0015 #include "TrackingTools/DetLayers/interface/ForwardDetLayer.h"
0016
0017 #include "RecoLocalTracker/SiStripClusterizer/interface/ClusterChargeCut.h"
0018
0019 #include "HitExtractorPIX.h"
0020 #include "HitExtractorSTRP.h"
0021
0022 #include <iostream>
0023 #include <sstream>
0024 #include <ostream>
0025 #include <fstream>
0026 #include <map>
0027
0028 using namespace ctfseeding;
0029 using namespace std;
0030
0031 SeedingLayerSetsBuilder::SeedingLayerId SeedingLayerSetsBuilder::nameToEnumId(const std::string& name) {
0032 GeomDetEnumerators::SubDetector subdet = GeomDetEnumerators::invalidDet;
0033 TrackerDetSide side = TrackerDetSide::Barrel;
0034 int idLayer = 0;
0035
0036 size_t index;
0037
0038
0039
0040 if ((index = name.find("BPix")) != string::npos) {
0041 subdet = GeomDetEnumerators::PixelBarrel;
0042 side = TrackerDetSide::Barrel;
0043 idLayer = atoi(name.substr(index + 4, 1).c_str());
0044 }
0045
0046
0047
0048 else if ((index = name.find("FPix")) != string::npos) {
0049 subdet = GeomDetEnumerators::PixelEndcap;
0050 idLayer = atoi(name.substr(index + 4).c_str());
0051 if (name.find("pos") != string::npos) {
0052 side = TrackerDetSide::PosEndcap;
0053 } else {
0054 side = TrackerDetSide::NegEndcap;
0055 }
0056 }
0057
0058
0059
0060 else if ((index = name.find("TIB")) != string::npos) {
0061 subdet = GeomDetEnumerators::TIB;
0062 side = TrackerDetSide::Barrel;
0063 idLayer = atoi(name.substr(index + 3, 1).c_str());
0064 }
0065
0066
0067
0068 else if ((index = name.find("TID")) != string::npos) {
0069 subdet = GeomDetEnumerators::TID;
0070 idLayer = atoi(name.substr(index + 3, 1).c_str());
0071 if (name.find("pos") != string::npos) {
0072 side = TrackerDetSide::PosEndcap;
0073 } else {
0074 side = TrackerDetSide::NegEndcap;
0075 }
0076 }
0077
0078
0079
0080 else if ((index = name.find("TOB")) != string::npos) {
0081 subdet = GeomDetEnumerators::TOB;
0082 side = TrackerDetSide::Barrel;
0083 idLayer = atoi(name.substr(index + 3, 1).c_str());
0084 }
0085
0086
0087
0088 else if ((index = name.find("TEC")) != string::npos) {
0089 subdet = GeomDetEnumerators::TEC;
0090 idLayer = atoi(name.substr(index + 3, 1).c_str());
0091 if (name.find("pos") != string::npos) {
0092 side = TrackerDetSide::PosEndcap;
0093 } else {
0094 side = TrackerDetSide::NegEndcap;
0095 }
0096 }
0097 return std::make_tuple(subdet, side, idLayer);
0098 }
0099
0100 SeedingLayerSetsBuilder::LayerSpec::LayerSpec(unsigned short index,
0101 const std::string& layerName,
0102 const edm::ParameterSet& cfgLayer,
0103 edm::ConsumesCollector& iC)
0104 : nameIndex(index),
0105 hitBuilder(cfgLayer.getParameter<string>("TTRHBuilder")),
0106 hitBuilderToken(iC.esConsumes(edm::ESInputTag("", hitBuilder))) {
0107 usePixelHitProducer = false;
0108 if (cfgLayer.exists("HitProducer")) {
0109 pixelHitProducer = cfgLayer.getParameter<string>("HitProducer");
0110 usePixelHitProducer = true;
0111 }
0112
0113 bool skipClusters = cfgLayer.exists("skipClusters");
0114 if (skipClusters) {
0115 LogDebug("SeedingLayerSetsBuilder") << layerName << " ready for skipping";
0116 } else {
0117 LogDebug("SeedingLayerSetsBuilder") << layerName << " not skipping ";
0118 }
0119
0120 auto subdetData = nameToEnumId(layerName);
0121 subdet = std::get<0>(subdetData);
0122 side = std::get<1>(subdetData);
0123 idLayer = std::get<2>(subdetData);
0124 if (subdet == GeomDetEnumerators::PixelBarrel || subdet == GeomDetEnumerators::PixelEndcap) {
0125 extractor = std::make_unique<HitExtractorPIX>(side, idLayer, pixelHitProducer, iC);
0126 } else if (subdet != GeomDetEnumerators::invalidDet) {
0127 auto extr = std::make_unique<HitExtractorSTRP>(subdet, side, idLayer, clusterChargeCut(cfgLayer), iC);
0128 if (cfgLayer.exists("matchedRecHits")) {
0129 extr->useMatchedHits(cfgLayer.getParameter<edm::InputTag>("matchedRecHits"), iC);
0130 }
0131 if (cfgLayer.exists("rphiRecHits")) {
0132 extr->useRPhiHits(cfgLayer.getParameter<edm::InputTag>("rphiRecHits"), iC);
0133 }
0134 if (cfgLayer.exists("stereoRecHits")) {
0135 extr->useStereoHits(cfgLayer.getParameter<edm::InputTag>("stereoRecHits"), iC);
0136 }
0137 if (cfgLayer.exists("vectorRecHits")) {
0138 extr->useVectorHits(cfgLayer.getParameter<edm::InputTag>("vectorRecHits"), iC);
0139 }
0140 if (cfgLayer.exists("useRingSlector") && cfgLayer.getParameter<bool>("useRingSlector")) {
0141 extr->useRingSelector(cfgLayer.getParameter<int>("minRing"), cfgLayer.getParameter<int>("maxRing"));
0142 }
0143 bool useSimpleRphiHitsCleaner =
0144 cfgLayer.exists("useSimpleRphiHitsCleaner") ? cfgLayer.getParameter<bool>("useSimpleRphiHitsCleaner") : true;
0145 extr->useSimpleRphiHitsCleaner(useSimpleRphiHitsCleaner);
0146
0147 double minAbsZ = cfgLayer.exists("MinAbsZ") ? cfgLayer.getParameter<double>("MinAbsZ") : 0.;
0148 if (minAbsZ > 0.) {
0149 extr->setMinAbsZ(minAbsZ);
0150 }
0151 if (skipClusters) {
0152 bool useProjection = cfgLayer.exists("useProjection") ? cfgLayer.getParameter<bool>("useProjection") : false;
0153 if (useProjection) {
0154 LogDebug("SeedingLayerSetsBuilder") << layerName << " will project partially masked matched rechit";
0155 } else {
0156 extr->setNoProjection();
0157 }
0158 }
0159 extractor = std::move(extr);
0160 }
0161 if (extractor && skipClusters) {
0162 extractor->useSkipClusters(cfgLayer.getParameter<edm::InputTag>("skipClusters"), iC);
0163 }
0164 }
0165
0166 std::string SeedingLayerSetsBuilder::LayerSpec::print(const std::vector<std::string>& names) const {
0167 std::ostringstream str;
0168 str << "Layer=" << names[nameIndex] << ", hitBldr: " << hitBuilder;
0169
0170 str << ", useRingSelector: ";
0171 HitExtractorSTRP* ext = nullptr;
0172 if ((ext = dynamic_cast<HitExtractorSTRP*>(extractor.get())) && ext->useRingSelector()) {
0173 auto minMaxRing = ext->getMinMaxRing();
0174 str << "true,"
0175 << " Rings: (" << std::get<0>(minMaxRing) << "," << std::get<1>(minMaxRing) << ")";
0176 } else
0177 str << "false";
0178
0179 return str.str();
0180 }
0181
0182 SeedingLayerSetsBuilder::SeedingLayerSetsBuilder(const edm::ParameterSet& cfg,
0183 edm::ConsumesCollector& iC,
0184 const edm::InputTag& fastsimHitTag)
0185 : SeedingLayerSetsBuilder(cfg, iC) {
0186 fastSimrecHitsToken_ = iC.consumes<FastTrackerRecHitCollection>(fastsimHitTag);
0187 trackerTopologyToken_ = iC.esConsumes();
0188 }
0189 SeedingLayerSetsBuilder::SeedingLayerSetsBuilder(const edm::ParameterSet& cfg, edm::ConsumesCollector&& iC)
0190 : SeedingLayerSetsBuilder(cfg, iC) {}
0191 SeedingLayerSetsBuilder::SeedingLayerSetsBuilder(const edm::ParameterSet& cfg, edm::ConsumesCollector& iC)
0192 : trackerToken_(iC.esConsumes()) {
0193 std::vector<std::string> namesPset = cfg.getParameter<std::vector<std::string> >("layerList");
0194 std::vector<std::vector<std::string> > layerNamesInSets = this->layerNamesInSets(namesPset);
0195
0196 typedef std::vector<std::string>::const_iterator IS;
0197 typedef std::vector<std::vector<std::string> >::const_iterator IT;
0198 std::ostringstream str;
0199
0200
0201
0202
0203
0204
0205
0206 if (layerNamesInSets.empty())
0207 theNumberOfLayersInSet = 0;
0208 else
0209 theNumberOfLayersInSet = layerNamesInSets[0].size();
0210
0211 for (IT it = layerNamesInSets.begin(); it != layerNamesInSets.end(); it++) {
0212 if (it->size() != theNumberOfLayersInSet)
0213 throw cms::Exception("Configuration")
0214 << "Assuming all SeedingLayerSets to have same number of layers. LayerSet " << (it - layerNamesInSets.begin())
0215 << " has " << it->size() << " while 0th has " << theNumberOfLayersInSet;
0216 for (const std::string& layerName : *it) {
0217 auto found = std::find(theLayerNames.begin(), theLayerNames.end(), layerName);
0218 unsigned short layerIndex = 0;
0219 if (found != theLayerNames.end()) {
0220 layerIndex = found - theLayerNames.begin();
0221 } else {
0222 if (std::numeric_limits<unsigned short>::max() == theLayerNames.size()) {
0223 throw cms::Exception("Assert")
0224 << "Too many layers in " << __FILE__ << ":" << __LINE__
0225 << ", we may have to enlarge the index type from unsigned short to unsigned int";
0226 }
0227
0228 layerIndex = theLayers.size();
0229 theLayers.emplace_back(theLayerNames.size(), layerName, layerConfig(layerName, cfg), iC);
0230 theLayerNames.push_back(layerName);
0231 }
0232 theLayerSetIndices.push_back(layerIndex);
0233 }
0234 }
0235 theLayerDets.resize(theLayers.size());
0236 theTTRHBuilders.resize(theLayers.size());
0237
0238
0239
0240
0241
0242
0243 }
0244
0245 SeedingLayerSetsBuilder::~SeedingLayerSetsBuilder() {}
0246
0247 void SeedingLayerSetsBuilder::fillDescriptions(edm::ParameterSetDescription& desc) {
0248 edm::ParameterSetDescription empty;
0249 empty.setAllowAnything();
0250
0251 desc.add<std::vector<std::string> >("layerList", {});
0252 desc.add<edm::ParameterSetDescription>("BPix", empty);
0253 desc.add<edm::ParameterSetDescription>("FPix", empty);
0254 desc.add<edm::ParameterSetDescription>("TIB", empty);
0255 desc.add<edm::ParameterSetDescription>("TID", empty);
0256 desc.add<edm::ParameterSetDescription>("TOB", empty);
0257 desc.add<edm::ParameterSetDescription>("TEC", empty);
0258 desc.add<edm::ParameterSetDescription>("MTIB", empty);
0259 desc.add<edm::ParameterSetDescription>("MTID", empty);
0260 desc.add<edm::ParameterSetDescription>("MTOB", empty);
0261 desc.add<edm::ParameterSetDescription>("MTEC", empty);
0262 }
0263
0264 edm::ParameterSet SeedingLayerSetsBuilder::layerConfig(const std::string& nameLayer,
0265 const edm::ParameterSet& cfg) const {
0266 edm::ParameterSet result;
0267
0268 for (string::size_type iEnd = nameLayer.size(); iEnd > 0; --iEnd) {
0269 string name = nameLayer.substr(0, iEnd);
0270 if (cfg.exists(name))
0271 return cfg.getParameter<edm::ParameterSet>(name);
0272 }
0273 edm::LogError("SeedingLayerSetsBuilder")
0274 << "configuration for layer: " << nameLayer << " not found, job will probably crash!";
0275 return result;
0276 }
0277
0278 vector<vector<string> > SeedingLayerSetsBuilder::layerNamesInSets(const vector<string>& namesPSet) {
0279 std::vector<std::vector<std::string> > result;
0280 for (std::vector<std::string>::const_iterator is = namesPSet.begin(); is < namesPSet.end(); ++is) {
0281 vector<std::string> layersInSet;
0282 string line = *is;
0283 string::size_type pos = 0;
0284 while (pos != string::npos) {
0285 pos = line.find('+');
0286 string layer = line.substr(0, pos);
0287 layersInSet.push_back(layer);
0288 line = line.substr(pos + 1, string::npos);
0289 }
0290 result.push_back(layersInSet);
0291 }
0292 return result;
0293 }
0294
0295 void SeedingLayerSetsBuilder::updateEventSetup(const edm::EventSetup& es) {
0296
0297
0298
0299 if (!(geometryWatcher_.check(es) || trhWatcher_.check(es)))
0300 return;
0301
0302 const GeometricSearchTracker& tracker = es.getData(trackerToken_);
0303
0304 const std::vector<BarrelDetLayer const*>& bpx = tracker.barrelLayers();
0305 const std::vector<BarrelDetLayer const*>& tib = tracker.tibLayers();
0306 const std::vector<BarrelDetLayer const*>& tob = tracker.tobLayers();
0307
0308 const std::vector<ForwardDetLayer const*>& fpx_pos = tracker.posForwardLayers();
0309 const std::vector<ForwardDetLayer const*>& tid_pos = tracker.posTidLayers();
0310 const std::vector<ForwardDetLayer const*>& tec_pos = tracker.posTecLayers();
0311
0312 const std::vector<ForwardDetLayer const*>& fpx_neg = tracker.negForwardLayers();
0313 const std::vector<ForwardDetLayer const*>& tid_neg = tracker.negTidLayers();
0314 const std::vector<ForwardDetLayer const*>& tec_neg = tracker.negTecLayers();
0315
0316 for (const auto& layer : theLayers) {
0317 const DetLayer* detLayer = nullptr;
0318 int index = layer.idLayer - 1;
0319
0320 if (layer.subdet == GeomDetEnumerators::PixelBarrel) {
0321 detLayer = bpx[index];
0322 } else if (layer.subdet == GeomDetEnumerators::PixelEndcap) {
0323 if (layer.side == TrackerDetSide::PosEndcap) {
0324 detLayer = fpx_pos[index];
0325 } else {
0326 detLayer = fpx_neg[index];
0327 }
0328 } else if (layer.subdet == GeomDetEnumerators::TIB) {
0329 detLayer = tib[index];
0330 } else if (layer.subdet == GeomDetEnumerators::TID) {
0331 if (layer.side == TrackerDetSide::PosEndcap) {
0332 detLayer = tid_pos[index];
0333 } else {
0334 detLayer = tid_neg[index];
0335 }
0336 } else if (layer.subdet == GeomDetEnumerators::TOB) {
0337 detLayer = tob[index];
0338 } else if (layer.subdet == GeomDetEnumerators::TEC) {
0339 if (layer.side == TrackerDetSide::PosEndcap) {
0340 detLayer = tec_pos[index];
0341 } else {
0342 detLayer = tec_neg[index];
0343 }
0344 } else {
0345 throw cms::Exception("Configuration") << "Did not find DetLayer for layer " << theLayerNames[layer.nameIndex];
0346 }
0347
0348 theLayerDets[layer.nameIndex] = detLayer;
0349 theTTRHBuilders[layer.nameIndex] = &es.getData(layer.hitBuilderToken);
0350 }
0351 }
0352
0353 std::vector<SeedingLayerSetsBuilder::SeedingLayerId> SeedingLayerSetsBuilder::layers() const {
0354 std::vector<SeedingLayerId> ret;
0355 ret.reserve(numberOfLayers());
0356 for (const auto& layer : theLayers) {
0357 ret.emplace_back(layer.subdet, layer.side, layer.idLayer);
0358 }
0359 return ret;
0360 }
0361
0362 std::unique_ptr<SeedingLayerSetsHits> SeedingLayerSetsBuilder::hits(const edm::Event& ev, const edm::EventSetup& es) {
0363 updateEventSetup(es);
0364
0365 auto ret = std::make_unique<SeedingLayerSetsHits>(
0366 theNumberOfLayersInSet, &theLayerSetIndices, &theLayerNames, &theLayerDets);
0367
0368 for (auto& layer : theLayers) {
0369 ret->addHits(
0370 layer.nameIndex,
0371 layer.extractor->hits((const TkTransientTrackingRecHitBuilder&)(*theTTRHBuilders[layer.nameIndex]), ev, es));
0372 }
0373 ret->shrink_to_fit();
0374 return ret;
0375 }
0376
0377 std::unique_ptr<SeedingLayerSetsHits> SeedingLayerSetsBuilder::makeSeedingLayerSetsHitsforFastSim(
0378 const edm::Event& ev, const edm::EventSetup& es) {
0379 updateEventSetup(es);
0380
0381 edm::Handle<FastTrackerRecHitCollection> fastSimrechits_;
0382 ev.getByToken(fastSimrecHitsToken_, fastSimrechits_);
0383 const TrackerTopology* const tTopo = &es.getData(trackerTopologyToken_);
0384 SeedingLayerSetsHits::OwnedHits layerhits_;
0385
0386 auto ret = std::make_unique<SeedingLayerSetsHits>(
0387 theNumberOfLayersInSet, &theLayerSetIndices, &theLayerNames, &theLayerDets);
0388
0389 for (auto& layer : theLayers) {
0390 layerhits_.clear();
0391 for (auto& rh : *fastSimrechits_) {
0392 GeomDetEnumerators::SubDetector subdet = GeomDetEnumerators::invalidDet;
0393 TrackerDetSide side = TrackerDetSide::Barrel;
0394 int idLayer = 0;
0395 if ((rh.det()->geographicalId()).subdetId() == PixelSubdetector::PixelBarrel) {
0396 subdet = GeomDetEnumerators::PixelBarrel;
0397 side = TrackerDetSide::Barrel;
0398 idLayer = tTopo->pxbLayer(rh.det()->geographicalId());
0399 } else if ((rh.det()->geographicalId()).subdetId() == PixelSubdetector::PixelEndcap) {
0400 subdet = GeomDetEnumerators::PixelEndcap;
0401 idLayer = tTopo->pxfDisk(rh.det()->geographicalId());
0402 if (tTopo->pxfSide(rh.det()->geographicalId()) == 1)
0403 side = TrackerDetSide::NegEndcap;
0404 else
0405 side = TrackerDetSide::PosEndcap;
0406 }
0407
0408 if (layer.subdet == subdet && layer.side == side && layer.idLayer == idLayer) {
0409 BaseTrackerRecHit const& b(rh);
0410 auto ptrHit = (BaseTrackerRecHit*)(b.clone());
0411 layerhits_.emplace_back(ptrHit);
0412 } else
0413 continue;
0414 }
0415 ret->addHits(layer.nameIndex, std::move(layerhits_));
0416 }
0417 ret->shrink_to_fit();
0418 return ret;
0419 }