MuonHits/src/GEMSimHitMatcher.cc

0001 #include "Validation/MuonHits/interface/GEMSimHitMatcher.h"
0002
0003 using namespace std;
0004
0005 GEMSimHitMatcher::GEMSimHitMatcher(const edm::ParameterSet& ps, edm::ConsumesCollector&& iC)
0006     : MuonSimHitMatcher(ps, std::move(iC)) {
0007   simHitPSet_ = ps.getParameterSet("gemSimHit");
0008   verbose_ = simHitPSet_.getParameter<int>("verbose");
0009   simMuOnly_ = simHitPSet_.getParameter<bool>("simMuOnly");
0010   discardEleHits_ = simHitPSet_.getParameter<bool>("discardEleHits");
0011
0012   simHitInput_ = iC.consumes<edm::PSimHitContainer>(simHitPSet_.getParameter<edm::InputTag>("inputTag"));
0013   geomToken_ = iC.esConsumes<GEMGeometry, MuonGeometryRecord>();
0014 }
0015
0016 /// initialize the event
0017 void GEMSimHitMatcher::init(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0018   geometry_ = &iSetup.getData(geomToken_);
0019   MuonSimHitMatcher::init(iEvent, iSetup);
0020 }
0021
0022 /// do the matching
0023 void GEMSimHitMatcher::match(const SimTrack& track, const SimVertex& vertex) {
0024   clear();
0025
0026   // instantiates the track ids and simhits
0027   MuonSimHitMatcher::match(track, vertex);
0028
0029   // hard cut on non-GEM muons
0030   if (std::abs(track.momentum().eta()) < 1.55)
0031     return;
0032
0033   matchSimHitsToSimTrack();
0034
0035   if (verbose_) {
0036     edm::LogInfo("GEMSimHitMatcher") << "nTrackIds " << track_ids_.size() << " nSelectedGEMSimHits " << hits_.size();
0037     edm::LogInfo("GEMSimHitMatcher") << "detids GEM " << detIds(0).size();
0038
0039     const auto& gem_ch_ids = detIds();
0040     for (const auto& id : gem_ch_ids) {
0041       const auto& gem_simhits = MuonSimHitMatcher::hitsInDetId(id);
0042       const auto& gem_simhits_gp = simHitsMeanPosition(gem_simhits);
0043       edm::LogInfo("GEMSimHitMatcher") << "gemchid " << GEMDetId(id) << ": nHits " << gem_simhits.size() << " phi "
0044                                        << gem_simhits_gp.phi() << " nCh " << chamber_to_hits_[id].size();
0045       const auto& strips = hitStripsInDetId(id);
0046       edm::LogInfo("GEMSimHitMatcher") << "nStrip " << strips.size();
0047       edm::LogInfo("GEMSimHitMatcher") << "strips : ";
0048       for (const auto& p : strips) {
0049         edm::LogInfo("GEMSimHitMatcher") << p;
0050       }
0051     }
0052     const auto& gem_sch_ids = superChamberIds();
0053     for (const auto& id : gem_sch_ids) {
0054       const auto& gem_simhits = hitsInSuperChamber(id);
0055       const auto& gem_simhits_gp = simHitsMeanPosition(gem_simhits);
0056       edm::LogInfo("GEMSimHitMatcher") << "gemschid " << GEMDetId(id) << ": " << nCoincidencePadsWithHits() << " | "
0057                                        << gem_simhits.size() << " " << gem_simhits_gp.phi() << " "
0058                                        << superchamber_to_hits_[id].size();
0059     }
0060   }
0061 }
0062
0063 void GEMSimHitMatcher::matchSimHitsToSimTrack() {
0064   for (const auto& track_id : track_ids_) {
0065     for (const auto& h : simHits_) {
0066       if (h.trackId() != track_id)
0067         continue;
0068
0069       const GEMDetId& p_id(h.detUnitId());
0070
0071       int pdgid = h.particleType();
0072
0073       // consider only the muon hits
0074       if (simMuOnly_ && std::abs(pdgid) != 13)
0075         continue;
0076
0077       // discard electron hits in the GEM chambers
0078       if (discardEleHits_ && std::abs(pdgid) == 11)
0079         continue;
0080
0081       detid_to_hits_[p_id.rawId()].push_back(h);
0082       hits_.push_back(h);
0083       chamber_to_hits_[p_id.chamberId().rawId()].push_back(h);
0084       superchamber_to_hits_[p_id.superChamberId().rawId()].push_back(h);
0085     }
0086   }
0087
0088   // find pads with hits
0089   const auto& detids = detIds();
0090   // find 2-layer coincidence pads with hits
0091   for (const auto& d : detids) {
0092     GEMDetId id(d);
0093     const auto& hits = hitsInDetId(d);
0094     const auto& roll = dynamic_cast<const GEMGeometry*>(geometry_)->etaPartition(id);
0095     // int max_npads = roll->npads();
0096     set<int> pads;
0097     for (const auto& h : hits) {
0098       const LocalPoint& lp = h.entryPoint();
0099       pads.insert(static_cast<int>(roll->padTopology().channel(lp)));
0100     }
0101     detids_to_pads_[d] = pads;
0102   }
0103
0104   // find 2-layer coincidence pads with hits
0105   for (const auto& d : detids) {
0106     GEMDetId id1(d);
0107     if (id1.layer() != 1)
0108       continue;
0109
0110     // find pads with hits in layer1
0111     const auto& hits1 = hitsInDetId(d);
0112     const auto& roll1 = dynamic_cast<const GEMGeometry*>(geometry_)->etaPartition(id1);
0113     set<int> pads1;
0114     set<int> pads2;
0115     set<int> copads;
0116
0117     for (const auto& h : hits1) {
0118       const LocalPoint& lp = h.entryPoint();
0119       pads1.insert(static_cast<int>(roll1->padTopology().channel(lp)));
0120       if (verbose_)
0121         edm::LogInfo("GEMSimHitMatcher") << "GEMHits detid1 " << id1 << " pad1 "
0122                                          << static_cast<int>(roll1->padTopology().channel(lp)) << std::endl;
0123     }
0124
0125     // find pads with hits in layer2
0126     for (const auto& d2 : detids) {
0127       // staggered geometry???? improve here !!
0128       GEMDetId id2(d2);
0129       // does layer 2 has simhits?
0130       if (id2.layer() != 2 or id2.region() != id1.region() or id2.ring() != id1.ring() or
0131           id2.station() != id1.station() or abs(id2.roll() - id1.roll()) > 1)
0132         continue;
0133       const auto& hits2 = hitsInDetId(id2());
0134       const auto& roll2 = dynamic_cast<const GEMGeometry*>(geometry_)->etaPartition(id2);
0135       for (const auto& h : hits2) {
0136         const LocalPoint& lp = h.entryPoint();
0137         pads2.insert(static_cast<int>(roll2->padTopology().channel(lp)));
0138         if (verbose_)
0139           edm::LogInfo("GEMSimHitMatcher") << "GEMHits detid2 " << id2 << " pad2 "
0140                                            << static_cast<int>(roll2->padTopology().channel(lp)) << std::endl;
0141       }
0142     }
0143
0144     for (const auto& pad1 : pads1) {
0145       for (const auto& pad2 : pads2) {
0146         if (abs(pad1 - pad2) <= 2) {
0147           if (copads.find(pad1) == copads.end())
0148             copads.insert(pad1);
0149           if (copads.find(pad2) == copads.end())
0150             copads.insert(pad2);
0151         }
0152       }
0153     }
0154
0155     if (copads.empty())
0156       continue;
0157
0158     // detids here is layer1 id
0159     detids_to_copads_[d] = copads;
0160   }
0161 }
0162
0163 std::set<unsigned int> GEMSimHitMatcher::detIds(int gem_type) const {
0164   std::set<unsigned int> result;
0165   for (const auto& p : detid_to_hits_) {
0166     const auto& id = p.first;
0167     if (gem_type > 0) {
0168       GEMDetId detId(id);
0169       if (MuonHitHelper::toGEMType(detId.station(), detId.ring()) != gem_type)
0170         continue;
0171     }
0172     result.insert(id);
0173   }
0174   return result;
0175 }
0176
0177 std::set<unsigned int> GEMSimHitMatcher::detIdsCoincidences() const {
0178   std::set<unsigned int> result;
0179   for (const auto& p : detids_to_copads_)
0180     result.insert(p.first);
0181   return result;
0182 }
0183
0184 std::set<unsigned int> GEMSimHitMatcher::chamberIds(int gem_type) const {
0185   std::set<unsigned int> result;
0186   for (const auto& p : chamber_to_hits_) {
0187     const auto& id = p.first;
0188     if (gem_type > 0) {
0189       GEMDetId detId(id);
0190       if (MuonHitHelper::toGEMType(detId.station(), detId.ring()) != gem_type)
0191         continue;
0192     }
0193     result.insert(id);
0194   }
0195   return result;
0196 }
0197
0198 std::set<unsigned int> GEMSimHitMatcher::superChamberIds() const {
0199   std::set<unsigned int> result;
0200   for (const auto& p : superchamber_to_hits_)
0201     result.insert(p.first);
0202   return result;
0203 }
0204
0205 std::set<unsigned int> GEMSimHitMatcher::superChamberIdsCoincidences() const {
0206   std::set<unsigned int> result;
0207   for (const auto& p : detids_to_copads_) {
0208     const GEMDetId& p_id(p.first);
0209     result.insert(p_id.superChamberId().rawId());
0210   }
0211   return result;
0212 }
0213
0214 const edm::PSimHitContainer& GEMSimHitMatcher::hitsInSuperChamber(unsigned int detid) const {
0215   if (MuonHitHelper::isGEM(detid)) {
0216     const GEMDetId id(detid);
0217     if (superchamber_to_hits_.find(id.chamberId().rawId()) == superchamber_to_hits_.end())
0218       return no_hits_;
0219     return superchamber_to_hits_.at(id.chamberId().rawId());
0220   }
0221
0222   return no_hits_;
0223 }
0224
0225 int GEMSimHitMatcher::nLayersWithHitsInSuperChamber(unsigned int detid) const {
0226   set<int> layers_with_hits;
0227   const auto& hits = hitsInSuperChamber(detid);
0228   for (const auto& h : hits) {
0229     const GEMDetId& idd(h.detUnitId());
0230     layers_with_hits.insert(idd.layer());
0231   }
0232   return layers_with_hits.size();
0233 }
0234
0235 bool GEMSimHitMatcher::hitStation(int st, int nlayers) const {
0236   int nst = 0;
0237   for (const auto& ddt : chamberIds()) {
0238     const GEMDetId id(ddt);
0239     if (id.station() != st)
0240       continue;
0241
0242     const int nl(nLayersWithHitsInSuperChamber(id.rawId()));
0243     if (nl < nlayers)
0244       continue;
0245     ++nst;
0246   }
0247   return nst;
0248 }
0249
0250 int GEMSimHitMatcher::nStations(int nlayers) const { return (hitStation(1, nlayers) + hitStation(2, nlayers)); }
0251
0252 float GEMSimHitMatcher::simHitsGEMCentralPosition(const edm::PSimHitContainer& sim_hits) const {
0253   if (sim_hits.empty())
0254     return -0.0;  // point "zero"
0255
0256   float central = -0.0;
0257   size_t n = 0;
0258   for (const auto& h : sim_hits) {
0259     LocalPoint lp(0., 0., 0.);  // local central
0260     GlobalPoint gp;
0261     if (MuonHitHelper::isGEM(h.detUnitId())) {
0262       gp = geometry_->idToDet(h.detUnitId())->surface().toGlobal(lp);
0263     }
0264     central = gp.perp();
0265     if (n >= 1)
0266       edm::LogWarning("GEMSimHitMatcher") << "warning! find more than one simhits in GEM chamber " << std::endl;
0267     ++n;
0268   }
0269
0270   return central;
0271 }
0272
0273 float GEMSimHitMatcher::simHitsMeanStrip(const edm::PSimHitContainer& sim_hits) const {
0274   if (sim_hits.empty())
0275     return -1.f;
0276
0277   float sums = 0.f;
0278   size_t n = 0;
0279   for (const auto& h : sim_hits) {
0280     const LocalPoint& lp = h.entryPoint();
0281     const auto& d = h.detUnitId();
0282     sums += dynamic_cast<const GEMGeometry*>(geometry_)->etaPartition(d)->strip(lp);
0283     ++n;
0284   }
0285   if (n == 0)
0286     return -1.f;
0287   return sums / n;
0288 }
0289
0290 std::set<int> GEMSimHitMatcher::hitStripsInDetId(unsigned int detid, int margin_n_strips) const {
0291   set<int> result;
0292   const auto& simhits = MuonSimHitMatcher::hitsInDetId(detid);
0293   GEMDetId id(detid);
0294   int max_nstrips = dynamic_cast<const GEMGeometry*>(geometry_)->etaPartition(id)->nstrips();
0295   for (const auto& h : simhits) {
0296     const LocalPoint& lp = h.entryPoint();
0297     int central_strip =
0298         static_cast<int>(dynamic_cast<const GEMGeometry*>(geometry_)->etaPartition(id)->topology().channel(lp));
0299     int smin = central_strip - margin_n_strips;
0300     smin = (smin > 0) ? smin : 0;
0301     int smax = central_strip + margin_n_strips;
0302     smax = (smax <= max_nstrips) ? smax : max_nstrips;
0303     for (int ss = smin; ss <= smax; ++ss)
0304       result.insert(ss);
0305   }
0306   return result;
0307 }
0308
0309 std::set<int> GEMSimHitMatcher::hitPadsInDetId(unsigned int detid) const {
0310   set<int> none;
0311   if (detids_to_pads_.find(detid) == detids_to_pads_.end())
0312     return none;
0313   return detids_to_pads_.at(detid);
0314 }
0315
0316 std::set<int> GEMSimHitMatcher::hitCoPadsInDetId(unsigned int detid) const {
0317   set<int> none;
0318   if (detids_to_copads_.find(detid) == detids_to_copads_.end())
0319     return none;
0320   return detids_to_copads_.at(detid);
0321 }
0322
0323 std::set<int> GEMSimHitMatcher::hitPartitions() const {
0324   std::set<int> result;
0325
0326   const auto& detids = detIds();
0327   for (const auto& id : detids) {
0328     GEMDetId idd(id);
0329     result.insert(idd.roll());
0330   }
0331   return result;
0332 }
0333
0334 int GEMSimHitMatcher::nPadsWithHits() const {
0335   int result = 0;
0336   const auto& pad_ids = detIds();
0337   for (const auto& id : pad_ids) {
0338     result += hitPadsInDetId(id).size();
0339   }
0340   return result;
0341 }
0342
0343 int GEMSimHitMatcher::nCoincidencePadsWithHits() const {
0344   int result = 0;
0345   const auto& copad_ids = detIdsCoincidences();
0346   for (const auto& id : copad_ids) {
0347     result += hitCoPadsInDetId(id).size();
0348   }
0349   return result;
0350 }
0351
0352 void GEMSimHitMatcher::clear() {
0353   MuonSimHitMatcher::clear();
0354
0355   superchamber_to_hits_.clear();
0356   detids_to_pads_.clear();
0357   detids_to_copads_.clear();
0358 }