File indexing completed on 2024-04-06 12:32:52
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
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
0024 void RPCSimHitMatcher::match(const SimTrack& track, const SimVertex& vertex) {
0025
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
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 }