MuonGEMRecHits/src/GEMRecHitMatcher.cc

0001 #include <memory>
0002
0003 #include "Validation/MuonGEMRecHits/interface/GEMRecHitMatcher.h"
0004 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0005 #include "Geometry/GEMGeometry/interface/GEMGeometry.h"
0006
0007 using namespace std;
0008
0009 GEMRecHitMatcher::GEMRecHitMatcher(const edm::ParameterSet& pset, edm::ConsumesCollector&& iC) {
0010   const auto& gemRecHit = pset.getParameterSet("gemRecHit");
0011   minBX_ = gemRecHit.getParameter<int>("minBX");
0012   maxBX_ = gemRecHit.getParameter<int>("maxBX");
0013   verbose_ = gemRecHit.getParameter<int>("verbose");
0014
0015   // make a new digi matcher
0016   gemDigiMatcher_ = std::make_shared<GEMDigiMatcher>(pset, std::move(iC));
0017
0018   gemRecHitToken_ = iC.consumes<GEMRecHitCollection>(gemRecHit.getParameter<edm::InputTag>("inputTag"));
0019   geomToken_ = iC.esConsumes<GEMGeometry, MuonGeometryRecord>();
0020 }
0021
0022 void GEMRecHitMatcher::init(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0023   gemDigiMatcher_->init(iEvent, iSetup);
0024
0025   iEvent.getByToken(gemRecHitToken_, gemRecHitH_);
0026
0027   gemGeometry_ = &iSetup.getData(geomToken_);
0028 }
0029
0030 /// do the matching
0031 void GEMRecHitMatcher::match(const SimTrack& t, const SimVertex& v) {
0032   // match digis first
0033   gemDigiMatcher_->match(t, v);
0034
0035   // get the rechit collection
0036   const GEMRecHitCollection& gemRecHits = *gemRecHitH_.product();
0037
0038   // now match the rechits
0039   matchRecHitsToSimTrack(gemRecHits);
0040 }
0041
0042 void GEMRecHitMatcher::matchRecHitsToSimTrack(const GEMRecHitCollection& rechits) {
0043   // get the matched ids with digis
0044   const auto& det_ids = gemDigiMatcher_->detIdsDigi();
0045
0046   // loop on those ids
0047   for (auto id : det_ids) {
0048     // now check the digis in this detid
0049     const auto& hit_digis = gemDigiMatcher_->stripNumbersInDetId(id);
0050     if (verbose()) {
0051       cout << "hit_digis_fat ";
0052       copy(hit_digis.begin(), hit_digis.end(), ostream_iterator<int>(cout, " "));
0053       cout << endl;
0054     }
0055
0056     GEMDetId p_id(id);
0057     const auto& rechits_in_det = rechits.get(p_id);
0058
0059     for (auto d = rechits_in_det.first; d != rechits_in_det.second; ++d) {
0060       if (verbose())
0061         edm::LogInfo("GEMDigiMatcher") << "recHit " << p_id << " " << *d << endl;
0062
0063       // check that the rechit is within BX range
0064       if (d->BunchX() < minBX_ || d->BunchX() > maxBX_)
0065         continue;
0066
0067       int firstStrip = d->firstClusterStrip();
0068       int cls = d->clusterSize();
0069       bool stripFound = false;
0070
0071       // check that it matches a strip that was hit by digis from our track
0072       for (int i = firstStrip; i < (firstStrip + cls); i++) {
0073         if (hit_digis.find(i) != hit_digis.end())
0074           stripFound = true;
0075       }
0076
0077       // this rechit did not correspond with any previously matched digi
0078       if (!stripFound)
0079         continue;
0080       if (verbose())
0081         edm::LogInfo("GEMDigiMatcher") << "oki" << endl;
0082
0083       recHits_.push_back(*d);
0084       detid_to_recHits_[id].push_back(*d);
0085       chamber_to_recHits_[p_id.chamberId().rawId()].push_back(*d);
0086       superchamber_to_recHits_[p_id.superChamberId().rawId()].push_back(*d);
0087     }
0088   }
0089 }
0090
0091 std::set<unsigned int> GEMRecHitMatcher::detIds() const {
0092   std::set<unsigned int> result;
0093   for (const auto& p : detid_to_recHits_)
0094     result.insert(p.first);
0095   return result;
0096 }
0097
0098 std::set<unsigned int> GEMRecHitMatcher::chamberIds() const {
0099   std::set<unsigned int> result;
0100   for (const auto& p : chamber_to_recHits_)
0101     result.insert(p.first);
0102   return result;
0103 }
0104
0105 std::set<unsigned int> GEMRecHitMatcher::superChamberIds() const {
0106   std::set<unsigned int> result;
0107   for (const auto& p : superchamber_to_recHits_)
0108     result.insert(p.first);
0109   return result;
0110 }
0111
0112 const GEMRecHitContainer& GEMRecHitMatcher::recHitsInDetId(unsigned int detid) const {
0113   if (detid_to_recHits_.find(detid) == detid_to_recHits_.end())
0114     return no_recHits_;
0115   return detid_to_recHits_.at(detid);
0116 }
0117
0118 const GEMRecHitContainer& GEMRecHitMatcher::recHitsInChamber(unsigned int detid) const {
0119   if (chamber_to_recHits_.find(detid) == chamber_to_recHits_.end())
0120     return no_recHits_;
0121   return chamber_to_recHits_.at(detid);
0122 }
0123
0124 const GEMRecHitContainer& GEMRecHitMatcher::recHitsInSuperChamber(unsigned int detid) const {
0125   if (superchamber_to_recHits_.find(detid) == superchamber_to_recHits_.end())
0126     return no_recHits_;
0127   return superchamber_to_recHits_.at(detid);
0128 }
0129
0130 int GEMRecHitMatcher::nLayersWithRecHitsInSuperChamber(unsigned int detid) const {
0131   set<int> layers;
0132   for (const auto& d : recHitsInSuperChamber(detid)) {
0133     layers.insert(d.gemId().layer());
0134   }
0135   return layers.size();
0136 }
0137
0138 std::set<int> GEMRecHitMatcher::stripNumbersInDetId(unsigned int detid) const {
0139   set<int> result;
0140   for (const auto& d : recHitsInDetId(detid)) {
0141     // loop on all strips hit in this rechit
0142     for (int iStrip = d.firstClusterStrip(); iStrip < d.firstClusterStrip() + d.clusterSize(); ++iStrip) {
0143       result.insert(iStrip);
0144     }
0145   }
0146   return result;
0147 }
0148
0149 std::set<int> GEMRecHitMatcher::partitionNumbers() const {
0150   std::set<int> result;
0151
0152   for (auto id : detIds()) {
0153     GEMDetId idd(id);
0154     result.insert(idd.roll());
0155   }
0156   return result;
0157 }
0158
0159 GlobalPoint GEMRecHitMatcher::recHitPosition(const GEMRecHit& rechit) const {
0160   const GEMDetId& idd = rechit.gemId();
0161   const LocalPoint& lp = rechit.localPosition();
0162   return gemGeometry_->idToDet(idd)->surface().toGlobal(lp);
0163 }
0164
0165 GlobalPoint GEMRecHitMatcher::recHitMeanPosition(const GEMRecHitContainer& rechit) const {
0166   GlobalPoint point_zero;
0167   if (rechit.empty())
0168     return point_zero;  // point "zero"
0169
0170   float sumx, sumy, sumz;
0171   sumx = sumy = sumz = 0.f;
0172   size_t n = 0;
0173   for (const auto& d : rechit) {
0174     GlobalPoint gp = recHitPosition(d);
0175     if (gp == point_zero)
0176       continue;
0177
0178     sumx += gp.x();
0179     sumy += gp.y();
0180     sumz += gp.z();
0181     ++n;
0182   }
0183   if (n == 0)
0184     return GlobalPoint();
0185   return GlobalPoint(sumx / n, sumy / n, sumz / n);
0186 }
0187
0188 bool GEMRecHitMatcher::recHitInContainer(const GEMRecHit& rh, const GEMRecHitContainer& c) const {
0189   bool isSame = false;
0190   for (const auto& thisRH : c)
0191     if (areGEMRecHitSame(thisRH, rh))
0192       isSame = true;
0193   return isSame;
0194 }
0195
0196 bool GEMRecHitMatcher::isGEMRecHitMatched(const GEMRecHit& thisRh) const {
0197   return recHitInContainer(thisRh, recHits());
0198 }
0199
0200 bool GEMRecHitMatcher::areGEMRecHitSame(const GEMRecHit& l, const GEMRecHit& r) const {
0201   return l.localPosition() == r.localPosition() and l.BunchX() == r.BunchX();
0202 }