File indexing completed on 2023-03-17 11:28:12
0001 #include "Validation/MuonHits/interface/ME0SimHitMatcher.h"
0002
0003 using namespace std;
0004
0005 ME0SimHitMatcher::ME0SimHitMatcher(const edm::ParameterSet& ps, edm::ConsumesCollector&& iC)
0006 : MuonSimHitMatcher(ps, std::move(iC)) {
0007 simHitPSet_ = ps.getParameterSet("me0SimHit");
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<ME0Geometry, MuonGeometryRecord>();
0014 }
0015
0016
0017 void ME0SimHitMatcher::init(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0018 geometry_ = &iSetup.getData(geomToken_);
0019 MuonSimHitMatcher::init(iEvent, iSetup);
0020 }
0021
0022
0023 void ME0SimHitMatcher::match(const SimTrack& track, const SimVertex& vertex) {
0024
0025 MuonSimHitMatcher::match(track, vertex);
0026
0027 matchSimHitsToSimTrack();
0028
0029 if (verbose_) {
0030 edm::LogInfo("ME0SimHitMatcher") << "nTrackIds " << track_ids_.size() << " nSelectedME0SimHits " << hits_.size();
0031 edm::LogInfo("ME0SimHitMatcher") << "detids ME0 " << detIds().size();
0032
0033 const auto& me0_ch_ids = detIds();
0034 for (const auto& id : me0_ch_ids) {
0035 const auto& me0_simhits = MuonSimHitMatcher::hitsInChamber(id);
0036 const auto& me0_simhits_gp = simHitsMeanPosition(me0_simhits);
0037 edm::LogInfo("ME0SimHitMatcher") << "me0chid " << ME0DetId(id) << ": nHits " << me0_simhits.size() << " phi "
0038 << me0_simhits_gp.phi() << " nCh " << chamber_to_hits_[id].size();
0039 const auto& strips = hitStripsInDetId(id);
0040 edm::LogInfo("ME0SimHitMatcher") << "nStrip " << strips.size();
0041 edm::LogInfo("ME0SimHitMatcher") << "strips : ";
0042 for (const auto& p : strips) {
0043 edm::LogInfo("ME0SimHitMatcher") << p;
0044 }
0045 }
0046 }
0047 }
0048
0049 void ME0SimHitMatcher::matchSimHitsToSimTrack() {
0050 for (const auto& track_id : track_ids_) {
0051 for (const auto& h : simHits_) {
0052 if (h.trackId() != track_id)
0053 continue;
0054 int pdgid = h.particleType();
0055 if (simMuOnly_ && std::abs(pdgid) != 13)
0056 continue;
0057
0058 if (discardEleHits_ && std::abs(pdgid) == 11)
0059 continue;
0060
0061 const ME0DetId& layer_id(h.detUnitId());
0062 detid_to_hits_[h.detUnitId()].push_back(h);
0063 hits_.push_back(h);
0064 chamber_to_hits_[layer_id.layerId().rawId()].push_back(h);
0065 superChamber_to_hits_[layer_id.chamberId().rawId()].push_back(h);
0066 }
0067 }
0068
0069
0070 const auto& detids = detIds();
0071 for (const auto& d : detids) {
0072 ME0DetId id(d);
0073 const auto& hits = hitsInDetId(d);
0074 const auto& roll = dynamic_cast<const ME0Geometry*>(geometry_)->etaPartition(id);
0075
0076 set<int> pads;
0077 for (const auto& h : hits) {
0078 const LocalPoint& lp = h.entryPoint();
0079 pads.insert(1 + static_cast<int>(roll->padTopology().channel(lp)));
0080 }
0081 detids_to_pads_[d] = pads;
0082 }
0083 }
0084
0085 std::set<unsigned int> ME0SimHitMatcher::detIds() const {
0086 std::set<unsigned int> result;
0087 for (const auto& p : detid_to_hits_)
0088 result.insert(p.first);
0089 return result;
0090 }
0091
0092 std::set<unsigned int> ME0SimHitMatcher::chamberIds() const {
0093 std::set<unsigned int> result;
0094 for (const auto& p : chamber_to_hits_)
0095 result.insert(p.first);
0096 return result;
0097 }
0098
0099 std::set<unsigned int> ME0SimHitMatcher::superChamberIds() const {
0100 std::set<unsigned int> result;
0101 for (const auto& p : superChamber_to_hits_)
0102 result.insert(p.first);
0103 return result;
0104 }
0105
0106 const edm::PSimHitContainer& ME0SimHitMatcher::hitsInSuperChamber(unsigned int detid) const {
0107 if (superChamber_to_hits_.find(detid) == superChamber_to_hits_.end())
0108 return no_hits_;
0109 return superChamber_to_hits_.at(detid);
0110 }
0111
0112 int ME0SimHitMatcher::nLayersWithHitsInSuperChamber(unsigned int detid) const {
0113 set<int> layers_with_hits;
0114 const auto& hits = hitsInSuperChamber(detid);
0115 for (const auto& h : hits) {
0116 const ME0DetId& idd(h.detUnitId());
0117 layers_with_hits.insert(idd.layer());
0118 }
0119 return layers_with_hits.size();
0120 }
0121
0122 std::set<unsigned int> ME0SimHitMatcher::superChamberIdsCoincidences(int min_n_layers) const {
0123 set<unsigned int> result;
0124
0125 for (const auto& p : superChamberIds()) {
0126 if (nLayersWithHitsInSuperChamber(p) >= min_n_layers) {
0127 result.insert(p);
0128 }
0129 }
0130 return result;
0131 }
0132
0133 int ME0SimHitMatcher::nCoincidenceChambers(int min_n_layers) const {
0134 return superChamberIdsCoincidences(min_n_layers).size();
0135 }
0136
0137 float ME0SimHitMatcher::simHitsMeanStrip(const edm::PSimHitContainer& sim_hits) const {
0138 if (sim_hits.empty())
0139 return -1.f;
0140
0141 float sums = 0.f;
0142 size_t n = 0;
0143 for (const auto& h : sim_hits) {
0144 const LocalPoint& lp = h.entryPoint();
0145 float s;
0146 const auto& d = h.detUnitId();
0147 s = dynamic_cast<const ME0Geometry*>(geometry_)->etaPartition(d)->strip(lp);
0148 sums += s;
0149 ++n;
0150 }
0151 if (n == 0)
0152 return -1.f;
0153 return sums / n;
0154 }
0155
0156 std::set<int> ME0SimHitMatcher::hitStripsInDetId(unsigned int detid, int margin_n_strips) const {
0157 set<int> result;
0158 const auto& simhits = hitsInDetId(detid);
0159 ME0DetId id(detid);
0160 int max_nstrips = dynamic_cast<const ME0Geometry*>(geometry_)->etaPartition(id)->nstrips();
0161 for (const auto& h : simhits) {
0162 const LocalPoint& lp = h.entryPoint();
0163 int central_strip =
0164 1 + static_cast<int>(dynamic_cast<const ME0Geometry*>(geometry_)->etaPartition(id)->topology().channel(lp));
0165 int smin = central_strip - margin_n_strips;
0166 smin = (smin > 0) ? smin : 1;
0167 int smax = central_strip + margin_n_strips;
0168 smax = (smax <= max_nstrips) ? smax : max_nstrips;
0169 for (int ss = smin; ss <= smax; ++ss)
0170 result.insert(ss);
0171 }
0172 return result;
0173 }
0174
0175 std::set<int> ME0SimHitMatcher::hitPadsInDetId(unsigned int detid) const {
0176 set<int> none;
0177 if (detids_to_pads_.find(detid) == detids_to_pads_.end())
0178 return none;
0179 return detids_to_pads_.at(detid);
0180 }
0181
0182 std::set<int> ME0SimHitMatcher::hitPartitions() const {
0183 std::set<int> result;
0184
0185 const auto& detids = detIds();
0186 for (const auto& id : detids) {
0187 ME0DetId idd(id);
0188 result.insert(idd.roll());
0189 }
0190 return result;
0191 }
0192
0193 int ME0SimHitMatcher::nPadsWithHits() const {
0194 int result = 0;
0195 const auto& pad_ids = detIds();
0196 for (const auto& id : pad_ids) {
0197 result += hitPadsInDetId(id).size();
0198 }
0199 return result;
0200 }