File indexing completed on 2025-01-22 07:34:36
0001 #include <memory>
0002
0003 #include "Validation/MuonGEMDigis/interface/GEMDigiMatcher.h"
0004
0005 using namespace std;
0006
0007 GEMDigiMatcher::GEMDigiMatcher(const edm::ParameterSet& pset, edm::ConsumesCollector&& iC) {
0008 const auto& gemSimLink = pset.getParameterSet("gemSimLink");
0009 simMuOnly_ = gemSimLink.getParameter<bool>("simMuOnly");
0010 discardEleHits_ = gemSimLink.getParameter<bool>("discardEleHits");
0011 verboseSimLink_ = gemSimLink.getParameter<int>("verbose");
0012
0013 const auto& gemDigi = pset.getParameterSet("gemStripDigi");
0014 minBXDigi_ = gemDigi.getParameter<int>("minBX");
0015 maxBXDigi_ = gemDigi.getParameter<int>("maxBX");
0016 matchDeltaStrip_ = gemDigi.getParameter<int>("matchDeltaStrip");
0017 verboseDigi_ = gemDigi.getParameter<int>("verbose");
0018 matchToSimLink_ = gemDigi.getParameter<bool>("matchToSimLink");
0019
0020 const auto& gemPad = pset.getParameterSet("gemPadDigi");
0021 minBXPad_ = gemPad.getParameter<int>("minBX");
0022 maxBXPad_ = gemPad.getParameter<int>("maxBX");
0023 verbosePad_ = gemPad.getParameter<int>("verbose");
0024
0025 const auto& gemCluster = pset.getParameterSet("gemPadCluster");
0026 minBXCluster_ = gemCluster.getParameter<int>("minBX");
0027 maxBXCluster_ = gemCluster.getParameter<int>("maxBX");
0028 verboseCluster_ = gemCluster.getParameter<int>("verbose");
0029
0030 const auto& gemCoPad = pset.getParameterSet("gemCoPadDigi");
0031 minBXCoPad_ = gemCoPad.getParameter<int>("minBX");
0032 maxBXCoPad_ = gemCoPad.getParameter<int>("maxBX");
0033 verboseCoPad_ = gemCoPad.getParameter<int>("verbose");
0034
0035
0036 muonSimHitMatcher_ = std::make_shared<GEMSimHitMatcher>(pset, std::move(iC));
0037
0038 if (matchToSimLink_)
0039 gemSimLinkToken_ =
0040 iC.consumes<edm::DetSetVector<GEMDigiSimLink>>(gemSimLink.getParameter<edm::InputTag>("inputTag"));
0041 gemDigiToken_ = iC.consumes<GEMDigiCollection>(gemDigi.getParameter<edm::InputTag>("inputTag"));
0042 gemPadToken_ = iC.consumes<GEMPadDigiCollection>(gemPad.getParameter<edm::InputTag>("inputTag"));
0043 gemClusterToken_ = iC.consumes<GEMPadDigiClusterCollection>(gemCluster.getParameter<edm::InputTag>("inputTag"));
0044 gemCoPadToken_ = iC.consumes<GEMCoPadDigiCollection>(gemCoPad.getParameter<edm::InputTag>("inputTag"));
0045
0046 geomToken_ = iC.esConsumes<GEMGeometry, MuonGeometryRecord>();
0047 }
0048
0049 void GEMDigiMatcher::init(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0050 muonSimHitMatcher_->init(iEvent, iSetup);
0051
0052 if (matchToSimLink_)
0053 iEvent.getByToken(gemSimLinkToken_, gemDigisSLH_);
0054 iEvent.getByToken(gemDigiToken_, gemDigisH_);
0055 iEvent.getByToken(gemPadToken_, gemPadsH_);
0056 iEvent.getByToken(gemClusterToken_, gemClustersH_);
0057 iEvent.getByToken(gemCoPadToken_, gemCoPadsH_);
0058
0059 const auto gemH = iSetup.getHandle(geomToken_);
0060 if (!gemH.isValid()) {
0061 gemGeometry_ = nullptr;
0062 edm::LogError("GEMDigiMatcher") << "Failed to initialize GEM geometry.";
0063 }
0064 gemGeometry_ = gemH.product();
0065 }
0066
0067
0068 void GEMDigiMatcher::match(const SimTrack& t, const SimVertex& v) {
0069
0070 muonSimHitMatcher_->match(t, v);
0071
0072
0073 const GEMDigiCollection& gemDigis = *gemDigisH_.product();
0074 const GEMPadDigiCollection& gemPads = *gemPadsH_.product();
0075 const GEMPadDigiClusterCollection& gemClusters = *gemClustersH_.product();
0076 const GEMCoPadDigiCollection& gemCoPads = *gemCoPadsH_.product();
0077
0078 clear();
0079
0080
0081 if (std::abs(t.momentum().eta()) < 1.55)
0082 return;
0083
0084
0085 if (matchToSimLink_) {
0086 const edm::DetSetVector<GEMDigiSimLink>& gemDigisSL = *gemDigisSLH_.product();
0087 matchDigisSLToSimTrack(gemDigisSL);
0088 }
0089 matchDigisToSimTrack(gemDigis);
0090 matchPadsToSimTrack(gemPads);
0091 matchClustersToSimTrack(gemClusters);
0092 matchCoPadsToSimTrack(gemCoPads);
0093 }
0094
0095 void GEMDigiMatcher::matchDigisSLToSimTrack(const edm::DetSetVector<GEMDigiSimLink>& digisSL) {
0096 if (verboseSimLink_)
0097 edm::LogInfo("GEMDigiMatcher") << "Matching simtrack to GEM simlinks" << endl;
0098
0099
0100 for (auto itsimlink = digisSL.begin(); itsimlink != digisSL.end(); itsimlink++) {
0101 GEMDetId p_id(itsimlink->id);
0102 for (auto sl = itsimlink->data.begin(); sl != itsimlink->data.end(); ++sl) {
0103
0104 const auto& detids(muonSimHitMatcher_->detIds());
0105 if (detids.find(p_id.rawId()) == detids.end())
0106 continue;
0107
0108
0109 if (muonSimHitMatcher_->hitsInDetId(p_id.rawId()).empty())
0110 continue;
0111
0112 if (verboseSimLink_)
0113 edm::LogInfo("GEMDigiMatcher") << "GEMDigiSimLink " << p_id << " " << sl->getStrip() << " " << sl->getBx()
0114 << " " << sl->getTrackId() << std::endl;
0115
0116
0117 if (simMuOnly_ && std::abs(sl->getParticleType()) != 13)
0118 continue;
0119
0120
0121 if (discardEleHits_ && std::abs(sl->getParticleType()) == 11)
0122 continue;
0123
0124
0125 for (const auto& simhit : muonSimHitMatcher_->hitsInDetId(p_id.rawId())) {
0126
0127 if (simhit.trackId() == sl->getTrackId() and simhit.particleType() == sl->getParticleType()) {
0128 detid_to_simLinks_[p_id.rawId()].push_back(*sl);
0129 if (verboseSimLink_)
0130 edm::LogInfo("GEMDigiMatcher") << "...was matched!" << endl;
0131 break;
0132 }
0133 }
0134 }
0135 }
0136 }
0137
0138 void GEMDigiMatcher::matchDigisToSimTrack(const GEMDigiCollection& digis) {
0139 if (verboseDigi_)
0140 edm::LogInfo("GEMDigiMatcher") << "Matching simtrack to GEM digis" << endl;
0141 for (auto id : muonSimHitMatcher_->detIds()) {
0142 GEMDetId p_id(id);
0143 const auto& hit_strips = muonSimHitMatcher_->hitStripsInDetId(id, matchDeltaStrip_);
0144 const auto& digis_in_det = digis.get(p_id);
0145
0146 for (auto d = digis_in_det.first; d != digis_in_det.second; ++d) {
0147 bool isMatched = false;
0148
0149
0150 if (d->bx() < minBXDigi_ || d->bx() > maxBXDigi_)
0151 continue;
0152
0153 if (verboseDigi_)
0154 edm::LogInfo("GEMDigiMatcher") << "GEMDigi " << p_id << " " << *d << endl;
0155
0156
0157 if (matchToSimLink_) {
0158
0159 for (const auto& sl : detid_to_simLinks_[p_id.rawId()]) {
0160 if (sl.getStrip() == d->strip() and sl.getBx() == d->bx()) {
0161 isMatched = true;
0162 break;
0163 }
0164 }
0165 }
0166
0167 else {
0168
0169 if (hit_strips.find(d->strip()) != hit_strips.end()) {
0170 isMatched = true;
0171 }
0172 }
0173 if (isMatched) {
0174 detid_to_digis_[p_id.rawId()].push_back(*d);
0175 chamber_to_digis_[p_id.chamberId().rawId()].push_back(*d);
0176 superchamber_to_digis_[p_id.superChamberId().rawId()].push_back(*d);
0177 if (verboseDigi_)
0178 edm::LogInfo("GEMDigiMatcher") << "...was matched!" << endl;
0179 }
0180 }
0181 }
0182 }
0183 void GEMDigiMatcher::matchPadsToSimTrack(const GEMPadDigiCollection& pads) {
0184 for (auto it = pads.begin(); it != pads.end(); ++it) {
0185 const GEMDetId& p_id = (*it).first;
0186 const auto& padsvec = (*it).second;
0187
0188 for (auto pad = padsvec.first; pad != padsvec.second; ++pad) {
0189
0190 if (pad->bx() < minBXPad_ || pad->bx() > maxBXPad_)
0191 continue;
0192
0193 if (verbosePad_)
0194 edm::LogInfo("GEMDigiMatcher") << "GEMPad " << p_id << " " << *pad << endl;
0195
0196 auto digivec = detid_to_digis_[p_id.rawId()];
0197
0198
0199
0200
0201
0202
0203
0204
0205
0206
0207
0208
0209 if (p_id.roll() % 2 == 1 && p_id.isGE21() && pad->nPartitions() == GEMPadDigi::GE21SplitStrip) {
0210
0211 GEMDetId p_id2(p_id.region(), p_id.ring(), p_id.station(), p_id.layer(), p_id.chamber(), p_id.roll() + 1);
0212
0213
0214 auto digivec2 = detid_to_digis_[p_id2.rawId()];
0215
0216
0217 digivec.insert(digivec.end(), digivec2.begin(), digivec2.end());
0218 }
0219 for (const auto& digi : digivec) {
0220
0221 const bool match8Partition(digi.strip() / 2 == pad->pad());
0222
0223
0224 const bool match16Partition(digi.strip() == pad->pad());
0225
0226
0227 const bool matchGE0(p_id.isME0() and match8Partition);
0228 const bool matchGE11(p_id.isGE11() and match8Partition);
0229 const bool matchGE21_8(p_id.isGE21() and pad->nPartitions() == GEMPadDigi::GE21 and match8Partition);
0230 const bool matchGE21_16(p_id.isGE21() and pad->nPartitions() == GEMPadDigi::GE21SplitStrip and
0231 match16Partition);
0232
0233
0234 if (matchGE0 or matchGE11 or matchGE21_8 or matchGE21_16) {
0235 detid_to_pads_[p_id.rawId()].push_back(*pad);
0236 chamber_to_pads_[p_id.chamberId().rawId()].push_back(*pad);
0237 superchamber_to_pads_[p_id.superChamberId().rawId()].push_back(*pad);
0238 if (verbosePad_)
0239 edm::LogInfo("GEMDigiMatcher") << "...was matched!" << endl;
0240 break;
0241 }
0242 }
0243 }
0244 }
0245 }
0246
0247 void GEMDigiMatcher::matchClustersToSimTrack(const GEMPadDigiClusterCollection& clusters) {
0248 for (auto it = clusters.begin(); it != clusters.end(); ++it) {
0249 const GEMDetId p_id = (*it).first;
0250 const auto clvec = (*it).second;
0251 for (auto cluster = clvec.first; cluster != clvec.second; ++cluster) {
0252 bool isMatched = false;
0253
0254
0255 if (cluster->bx() < minBXCluster_ || cluster->bx() > maxBXCluster_)
0256 continue;
0257
0258 if (verboseCluster_)
0259 edm::LogInfo("GEMDigiMatcher") << "GEMCluster " << p_id << " " << *cluster << endl;
0260
0261
0262 for (const auto& p : cluster->pads()) {
0263 for (const auto& pad : detid_to_pads_[p_id.rawId()]) {
0264 if (pad.pad() == p) {
0265 isMatched = true;
0266 }
0267 }
0268 }
0269 if (isMatched) {
0270 detid_to_clusters_[p_id.rawId()].push_back(*cluster);
0271 chamber_to_clusters_[p_id.chamberId().rawId()].push_back(*cluster);
0272 superchamber_to_clusters_[p_id.superChamberId().rawId()].push_back(*cluster);
0273 if (verboseCluster_)
0274 edm::LogInfo("GEMDigiMatcher") << "...was matched!" << endl;
0275 }
0276 }
0277 }
0278 }
0279
0280 void GEMDigiMatcher::matchCoPadsToSimTrack(const GEMCoPadDigiCollection& co_pads) {
0281
0282 for (auto d : superChamberIdsPad()) {
0283 GEMDetId id(d);
0284
0285 const auto& co_pads_in_det = co_pads.get(id);
0286 for (auto copad = co_pads_in_det.first; copad != co_pads_in_det.second; ++copad) {
0287
0288 if (copad->bx(1) < minBXCoPad_ || copad->bx(1) > maxBXCoPad_)
0289 continue;
0290
0291 if (verboseCoPad_)
0292 edm::LogInfo("GEMDigiMatcher") << "GEMCoPadDigi: " << id << " " << *copad << endl;
0293
0294 bool isMatchedL1 = false;
0295 bool isMatchedL2 = false;
0296 GEMDetId gemL1_id(id.region(), 1, id.station(), 1, id.chamber(), copad->roll());
0297 GEMDetId gemL2_id(id.region(), 1, id.station(), 2, id.chamber(), 0);
0298
0299
0300 for (const auto& p : padsInDetId(gemL1_id.rawId())) {
0301 if (p == copad->first()) {
0302 isMatchedL1 = true;
0303 }
0304 }
0305
0306
0307 for (const auto& p : padsInChamber(gemL2_id.rawId())) {
0308 if (p == copad->second()) {
0309 isMatchedL2 = true;
0310 }
0311 }
0312 if (isMatchedL1 and isMatchedL2) {
0313 superchamber_to_copads_[id.rawId()].push_back(*copad);
0314 if (verboseCoPad_)
0315 edm::LogInfo("GEMDigiMatcher") << "...was matched! " << endl;
0316 }
0317 }
0318 }
0319 }
0320
0321 std::set<unsigned int> GEMDigiMatcher::detIdsSimLink(int gem_type) const {
0322 return selectDetIds(detid_to_simLinks_, gem_type);
0323 }
0324
0325 std::set<unsigned int> GEMDigiMatcher::detIdsDigi(int gem_type) const {
0326 return selectDetIds(detid_to_digis_, gem_type);
0327 }
0328
0329 std::set<unsigned int> GEMDigiMatcher::detIdsPad(int gem_type) const { return selectDetIds(detid_to_pads_, gem_type); }
0330
0331 std::set<unsigned int> GEMDigiMatcher::detIdsCluster(int gem_type) const {
0332 return selectDetIds(detid_to_clusters_, gem_type);
0333 }
0334
0335 std::set<unsigned int> GEMDigiMatcher::chamberIdsDigi(int gem_type) const {
0336 return selectDetIds(chamber_to_digis_, gem_type);
0337 }
0338
0339 std::set<unsigned int> GEMDigiMatcher::chamberIdsPad(int gem_type) const {
0340 return selectDetIds(chamber_to_pads_, gem_type);
0341 }
0342
0343 std::set<unsigned int> GEMDigiMatcher::chamberIdsCluster(int gem_type) const {
0344 return selectDetIds(chamber_to_clusters_, gem_type);
0345 }
0346
0347 std::set<unsigned int> GEMDigiMatcher::superChamberIdsDigi(int gem_type) const {
0348 return selectDetIds(superchamber_to_digis_, gem_type);
0349 }
0350
0351 std::set<unsigned int> GEMDigiMatcher::superChamberIdsPad(int gem_type) const {
0352 return selectDetIds(superchamber_to_pads_, gem_type);
0353 }
0354
0355 std::set<unsigned int> GEMDigiMatcher::superChamberIdsCluster(int gem_type) const {
0356 return selectDetIds(superchamber_to_clusters_, gem_type);
0357 }
0358
0359 std::set<unsigned int> GEMDigiMatcher::superChamberIdsCoPad(int gem_type) const {
0360 return selectDetIds(superchamber_to_copads_, gem_type);
0361 }
0362
0363 const GEMDigiContainer& GEMDigiMatcher::digisInDetId(unsigned int detid) const {
0364 if (detid_to_digis_.find(detid) == detid_to_digis_.end())
0365 return no_gem_digis_;
0366 return detid_to_digis_.at(detid);
0367 }
0368
0369 const GEMDigiContainer& GEMDigiMatcher::digisInChamber(unsigned int detid) const {
0370 if (chamber_to_digis_.find(detid) == chamber_to_digis_.end())
0371 return no_gem_digis_;
0372 return chamber_to_digis_.at(detid);
0373 }
0374
0375 const GEMDigiContainer& GEMDigiMatcher::digisInSuperChamber(unsigned int detid) const {
0376 if (superchamber_to_digis_.find(detid) == superchamber_to_digis_.end())
0377 return no_gem_digis_;
0378 return superchamber_to_digis_.at(detid);
0379 }
0380
0381 const GEMPadDigiContainer& GEMDigiMatcher::padsInDetId(unsigned int detid) const {
0382 if (detid_to_pads_.find(detid) == detid_to_pads_.end())
0383 return no_gem_pads_;
0384 return detid_to_pads_.at(detid);
0385 }
0386
0387 const GEMPadDigiContainer& GEMDigiMatcher::padsInChamber(unsigned int detid) const {
0388 if (chamber_to_pads_.find(detid) == chamber_to_pads_.end())
0389 return no_gem_pads_;
0390 return chamber_to_pads_.at(detid);
0391 }
0392
0393 const GEMPadDigiContainer& GEMDigiMatcher::padsInSuperChamber(unsigned int detid) const {
0394 if (superchamber_to_pads_.find(detid) == superchamber_to_pads_.end())
0395 return no_gem_pads_;
0396 return superchamber_to_pads_.at(detid);
0397 }
0398
0399 const GEMPadDigiClusterContainer& GEMDigiMatcher::clustersInDetId(unsigned int detid) const {
0400 if (detid_to_clusters_.find(detid) == detid_to_clusters_.end())
0401 return no_gem_clusters_;
0402 return detid_to_clusters_.at(detid);
0403 }
0404
0405 const GEMPadDigiClusterContainer& GEMDigiMatcher::clustersInChamber(unsigned int detid) const {
0406 if (chamber_to_clusters_.find(detid) == chamber_to_clusters_.end())
0407 return no_gem_clusters_;
0408 return chamber_to_clusters_.at(detid);
0409 }
0410
0411 const GEMPadDigiClusterContainer& GEMDigiMatcher::clustersInSuperChamber(unsigned int detid) const {
0412 if (superchamber_to_clusters_.find(detid) == superchamber_to_clusters_.end())
0413 return no_gem_clusters_;
0414 return superchamber_to_clusters_.at(detid);
0415 }
0416
0417 const GEMCoPadDigiContainer& GEMDigiMatcher::coPadsInSuperChamber(unsigned int detid) const {
0418 if (superchamber_to_copads_.find(detid) == superchamber_to_copads_.end())
0419 return no_gem_copads_;
0420 return superchamber_to_copads_.at(detid);
0421 }
0422
0423 int GEMDigiMatcher::nLayersWithDigisInSuperChamber(unsigned int detid) const {
0424 set<int> layers;
0425 GEMDetId sch_id(detid);
0426 for (int iLayer = 1; iLayer <= 2; iLayer++) {
0427 GEMDetId ch_id(sch_id.region(), sch_id.ring(), sch_id.station(), iLayer, sch_id.chamber(), 0);
0428
0429 const auto& digis = digisInChamber(ch_id.rawId());
0430
0431 if (!digis.empty()) {
0432 layers.insert(iLayer);
0433 }
0434 }
0435 return layers.size();
0436 }
0437
0438 int GEMDigiMatcher::nLayersWithPadsInSuperChamber(unsigned int detid) const {
0439 set<int> layers;
0440 GEMDetId sch_id(detid);
0441 for (int iLayer = 1; iLayer <= 2; iLayer++) {
0442 GEMDetId ch_id(sch_id.region(), sch_id.ring(), sch_id.station(), iLayer, sch_id.chamber(), 0);
0443
0444 const auto& pads = padsInChamber(ch_id.rawId());
0445
0446 if (!pads.empty()) {
0447 layers.insert(iLayer);
0448 }
0449 }
0450 return layers.size();
0451 }
0452
0453 int GEMDigiMatcher::nLayersWithClustersInSuperChamber(unsigned int detid) const {
0454 set<int> layers;
0455 GEMDetId sch_id(detid);
0456 for (int iLayer = 1; iLayer <= 2; iLayer++) {
0457 GEMDetId ch_id(sch_id.region(), sch_id.ring(), sch_id.station(), iLayer, sch_id.chamber(), 0);
0458
0459 const auto& clusters = clustersInChamber(ch_id.rawId());
0460
0461 if (!clusters.empty()) {
0462 layers.insert(iLayer);
0463 }
0464 }
0465 return layers.size();
0466 }
0467
0468 int GEMDigiMatcher::nPads() const {
0469 int n = 0;
0470 const auto& ids = superChamberIdsPad();
0471 for (const auto& id : ids) {
0472 n += padsInSuperChamber(id).size();
0473 }
0474 return n;
0475 }
0476
0477 int GEMDigiMatcher::nCoPads() const {
0478 int n = 0;
0479 const auto& ids = superChamberIdsCoPad();
0480 for (const auto& id : ids) {
0481 n += coPadsInSuperChamber(id).size();
0482 }
0483 return n;
0484 }
0485
0486 std::set<int> GEMDigiMatcher::stripNumbersInDetId(unsigned int detid) const {
0487 set<int> result;
0488 const auto& digis = digisInDetId(detid);
0489 for (const auto& d : digis) {
0490 result.insert(d.strip());
0491 }
0492 return result;
0493 }
0494
0495 std::set<int> GEMDigiMatcher::padNumbersInDetId(unsigned int detid) const {
0496 set<int> result;
0497 const auto& digis = padsInDetId(detid);
0498 for (const auto& d : digis) {
0499 result.insert(d.pad());
0500 }
0501 return result;
0502 }
0503
0504 std::set<int> GEMDigiMatcher::partitionNumbers() const {
0505 std::set<int> result;
0506
0507 const auto& detids = detIdsDigi();
0508 for (const auto& id : detids) {
0509 const GEMDetId& idd(id);
0510 result.insert(idd.roll());
0511 }
0512 return result;
0513 }
0514
0515 std::set<int> GEMDigiMatcher::partitionNumbersWithCoPads() const {
0516 std::set<int> result;
0517
0518 const auto& detids = superChamberIdsCoPad();
0519 for (const auto& id : detids) {
0520 const GEMDetId& idd(id);
0521 result.insert(idd.roll());
0522 }
0523 return result;
0524 }
0525
0526 GlobalPoint GEMDigiMatcher::getGlobalPointDigi(unsigned int rawId, const GEMDigi& d) const {
0527 GEMDetId gem_id(rawId);
0528 const LocalPoint& gem_lp = gemGeometry_->etaPartition(gem_id)->centreOfStrip(d.strip());
0529 const GlobalPoint& gem_gp = gemGeometry_->idToDet(gem_id)->surface().toGlobal(gem_lp);
0530 return gem_gp;
0531 }
0532
0533 GlobalPoint GEMDigiMatcher::getGlobalPointPad(unsigned int rawId, const GEMPadDigi& tp) const {
0534 GEMDetId gem_id(rawId);
0535 const LocalPoint& gem_lp = gemGeometry_->etaPartition(gem_id)->centreOfPad(tp.pad());
0536 const GlobalPoint& gem_gp = gemGeometry_->idToDet(gem_id)->surface().toGlobal(gem_lp);
0537 return gem_gp;
0538 }
0539
0540 void GEMDigiMatcher::clear() {
0541 detid_to_simLinks_.clear();
0542
0543 detid_to_digis_.clear();
0544 chamber_to_digis_.clear();
0545 superchamber_to_digis_.clear();
0546
0547 detid_to_pads_.clear();
0548 chamber_to_pads_.clear();
0549 superchamber_to_pads_.clear();
0550
0551 detid_to_clusters_.clear();
0552 chamber_to_clusters_.clear();
0553 superchamber_to_clusters_.clear();
0554
0555 superchamber_to_copads_.clear();
0556 }