MuonHits/src/RPCSimHitMatcher.cc

0001 #include "Validation/MuonHits/interface/RPCSimHitMatcher.h"
0002 #include "Geometry/CommonTopologies/interface/TrapezoidalStripTopology.h"
0003
0004 using namespace std;
0005
0006 RPCSimHitMatcher::RPCSimHitMatcher(const edm::ParameterSet& ps, edm::ConsumesCollector&& iC)
0007     : MuonSimHitMatcher(ps, std::move(iC)) {
0008   simHitPSet_ = ps.getParameterSet("rpcSimHit");
0009   verbose_ = simHitPSet_.getParameter<int>("verbose");
0010   simMuOnly_ = simHitPSet_.getParameter<bool>("simMuOnly");
0011   discardEleHits_ = simHitPSet_.getParameter<bool>("discardEleHits");
0012
0013   simHitInput_ = iC.consumes<edm::PSimHitContainer>(simHitPSet_.getParameter<edm::InputTag>("inputTag"));
0014   geomToken_ = iC.esConsumes<RPCGeometry, MuonGeometryRecord>();
0015 }
0016
0017 /// initialize the event
0018 void RPCSimHitMatcher::init(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0019   geometry_ = &iSetup.getData(geomToken_);
0020   MuonSimHitMatcher::init(iEvent, iSetup);
0021 }
0022
0023 /// do the matching
0024 void RPCSimHitMatcher::match(const SimTrack& track, const SimVertex& vertex) {
0025   // instantiates the track ids and simhits
0026   MuonSimHitMatcher::match(track, vertex);
0027
0028   matchSimHitsToSimTrack();
0029
0030   if (verbose_) {
0031     edm::LogInfo("RPCSimHitMatcher") << "nSimHits " << simHits_.size() << " nTrackIds " << track_ids_.size();
0032     edm::LogInfo("RPCSimHitMatcher") << "detids RPC " << detIds().size();
0033
0034     const auto& ch_ids = chamberIds();
0035     for (const auto& id : ch_ids) {
0036       const auto& simhits = MuonSimHitMatcher::hitsInChamber(id);
0037       const auto& simhits_gp = simHitsMeanPosition(simhits);
0038       edm::LogInfo("RPCSimHitMatcher") << "RPCDetId " << RPCDetId(id) << ": nHits " << simhits.size() << " eta "
0039                                        << simhits_gp.eta() << " phi " << simhits_gp.phi() << " nCh "
0040                                        << chamber_to_hits_[id].size();
0041       const auto& strips = hitStripsInDetId(id);
0042       edm::LogInfo("RPCSimHitMatcher") << "nStrips " << strips.size();
0043       edm::LogInfo("RPCSimHitMatcher") << "strips : ";
0044       for (const auto& p : strips) {
0045         edm::LogInfo("RPCSimHitMatcher") << p;
0046       }
0047     }
0048   }
0049 }
0050
0051 void RPCSimHitMatcher::matchSimHitsToSimTrack() {
0052   for (const auto& track_id : track_ids_) {
0053     for (const auto& h : simHits_) {
0054       if (h.trackId() != track_id)
0055         continue;
0056       int pdgid = h.particleType();
0057       if (simMuOnly_ && std::abs(pdgid) != 13)
0058         continue;
0059       // discard electron hits in the RPC chambers
0060       if (discardEleHits_ && std::abs(pdgid) == 11)
0061         continue;
0062
0063       const RPCDetId& layer_id(h.detUnitId());
0064       detid_to_hits_[h.detUnitId()].push_back(h);
0065       hits_.push_back(h);
0066       chamber_to_hits_[layer_id.chamberId().rawId()].push_back(h);
0067     }
0068   }
0069 }
0070
0071 std::set<unsigned int> RPCSimHitMatcher::detIds(int type) const {
0072   std::set<unsigned int> result;
0073   for (const auto& p : detid_to_hits_) {
0074     const auto& id = p.first;
0075     if (type > 0) {
0076       RPCDetId detId(id);
0077       if (MuonHitHelper::toRPCType(detId.region(), detId.station(), detId.ring()) != type)
0078         continue;
0079     }
0080     result.insert(id);
0081   }
0082   return result;
0083 }
0084
0085 std::set<unsigned int> RPCSimHitMatcher::chamberIds(int type) const {
0086   std::set<unsigned int> result;
0087   for (const auto& p : chamber_to_hits_) {
0088     const auto& id = p.first;
0089     if (type > 0) {
0090       RPCDetId detId(id);
0091       if (MuonHitHelper::toRPCType(detId.region(), detId.station(), detId.ring()) != type)
0092         continue;
0093     }
0094     result.insert(id);
0095   }
0096   return result;
0097 }
0098
0099 bool RPCSimHitMatcher::hitStation(int st) const {
0100   int nst = 0;
0101   for (const auto& ddt : chamberIds(0)) {
0102     const RPCDetId id(ddt);
0103     if (id.station() != st)
0104       continue;
0105     ++nst;
0106   }
0107   return nst;
0108 }
0109
0110 int RPCSimHitMatcher::nStations() const { return (hitStation(1) + hitStation(2) + hitStation(3) + hitStation(4)); }
0111
0112 float RPCSimHitMatcher::simHitsMeanStrip(const edm::PSimHitContainer& sim_hits) const {
0113   if (sim_hits.empty())
0114     return -1.f;
0115
0116   float sums = 0.f;
0117   size_t n = 0;
0118   for (const auto& h : sim_hits) {
0119     const LocalPoint& lp = h.entryPoint();
0120     const auto& d = h.detUnitId();
0121     sums += dynamic_cast<const RPCGeometry*>(geometry_)->roll(d)->strip(lp);
0122     ++n;
0123   }
0124   if (n == 0)
0125     return -1.f;
0126   return sums / n;
0127 }
0128
0129 std::set<int> RPCSimHitMatcher::hitStripsInDetId(unsigned int detid, int margin_n_strips) const {
0130   set<int> result;
0131   RPCDetId id(detid);
0132   for (const auto& roll : dynamic_cast<const RPCGeometry*>(geometry_)->chamber(id)->rolls()) {
0133     int max_nstrips = roll->nstrips();
0134     for (const auto& h : MuonSimHitMatcher::hitsInDetId(roll->id().rawId())) {
0135       const LocalPoint& lp = h.entryPoint();
0136       int central_strip = static_cast<int>(roll->topology().channel(lp));
0137       int smin = central_strip - margin_n_strips;
0138       smin = (smin > 0) ? smin : 1;
0139       int smax = central_strip + margin_n_strips;
0140       smax = (smax <= max_nstrips) ? smax : max_nstrips;
0141       for (int ss = smin; ss <= smax; ++ss)
0142         result.insert(ss);
0143     }
0144   }
0145   return result;
0146 }