MuonGEMRecHits/src/GEMRecHitMatcher.cc

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