File indexing completed on 2024-04-06 12:31:57
0001 #include <memory>
0002
0003 #include "Validation/CSCRecHits/interface/CSCRecHitMatcher.h"
0004
0005 using namespace std;
0006
0007 CSCRecHitMatcher::CSCRecHitMatcher(const edm::ParameterSet& pset, edm::ConsumesCollector&& iC) {
0008 const auto& cscRecHit2D = pset.getParameter<edm::ParameterSet>("cscRecHit");
0009 maxBXCSCRecHit2D_ = cscRecHit2D.getParameter<int>("maxBX");
0010 minBXCSCRecHit2D_ = cscRecHit2D.getParameter<int>("minBX");
0011 verboseCSCRecHit2D_ = cscRecHit2D.getParameter<int>("verbose");
0012
0013 const auto& cscSegment = pset.getParameter<edm::ParameterSet>("cscSegment");
0014 maxBXCSCSegment_ = cscSegment.getParameter<int>("maxBX");
0015 minBXCSCSegment_ = cscSegment.getParameter<int>("minBX");
0016 verboseCSCSegment_ = cscSegment.getParameter<int>("verbose");
0017
0018
0019 cscDigiMatcher_ = std::make_unique<CSCDigiMatcher>(pset, std::move(iC));
0020
0021 cscRecHit2DToken_ = iC.consumes<CSCRecHit2DCollection>(cscRecHit2D.getParameter<edm::InputTag>("inputTag"));
0022 cscSegmentToken_ = iC.consumes<CSCSegmentCollection>(cscSegment.getParameter<edm::InputTag>("inputTag"));
0023
0024 cscGeomToken_ = iC.esConsumes<CSCGeometry, MuonGeometryRecord>();
0025 }
0026
0027 void CSCRecHitMatcher::init(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0028 cscDigiMatcher_->init(iEvent, iSetup);
0029
0030 iEvent.getByToken(cscRecHit2DToken_, cscRecHit2DH_);
0031 iEvent.getByToken(cscSegmentToken_, cscSegmentH_);
0032
0033 const auto& csc_geom = iSetup.getHandle(cscGeomToken_);
0034 if (csc_geom.isValid()) {
0035 cscGeometry_ = csc_geom.product();
0036 } else {
0037 edm::LogWarning("CSCSimHitMatcher") << "+++ Info: CSC geometry is unavailable. +++\n";
0038 }
0039 }
0040
0041
0042 void CSCRecHitMatcher::match(const SimTrack& t, const SimVertex& v) {
0043
0044 cscDigiMatcher_->match(t, v);
0045
0046
0047 const CSCRecHit2DCollection& cscRecHit2Ds = *cscRecHit2DH_.product();
0048 const CSCSegmentCollection& cscSegments = *cscSegmentH_.product();
0049
0050
0051 matchCSCRecHit2DsToSimTrack(cscRecHit2Ds);
0052 matchCSCSegmentsToSimTrack(cscSegments);
0053 }
0054
0055 void CSCRecHitMatcher::matchCSCRecHit2DsToSimTrack(const CSCRecHit2DCollection& rechits) {
0056 if (verboseCSCRecHit2D_)
0057 edm::LogInfo("CSCRecHitMatcher") << "Matching simtrack to CSC rechits";
0058
0059
0060 const auto& strip_ids = cscDigiMatcher_->detIdsStrip();
0061 const auto& wire_ids = cscDigiMatcher_->detIdsWire();
0062
0063
0064 std::set<unsigned int> layer_ids;
0065 layer_ids.insert(strip_ids.begin(), strip_ids.end());
0066 layer_ids.insert(wire_ids.begin(), wire_ids.end());
0067
0068 if (verboseCSCRecHit2D_)
0069 edm::LogInfo("CSCRecHitMatcher") << "Number of matched csc layer_ids " << layer_ids.size();
0070
0071 for (const auto& id : layer_ids) {
0072 CSCDetId p_id(id);
0073
0074
0075 const auto& hit_wg(cscDigiMatcher_->wiregroupsInDetId(id));
0076 if (verboseCSCRecHit2D_) {
0077 edm::LogInfo("CSCRecHitMatcher") << "hit wg csc from simhit" << endl;
0078 for (const auto& p : hit_wg) {
0079 edm::LogInfo("CSCRecHitMatcher") << p;
0080 }
0081 }
0082
0083
0084 const auto& hit_strips(cscDigiMatcher_->stripsInDetId(id));
0085 if (verboseCSCRecHit2D_) {
0086 edm::LogInfo("CSCRecHitMatcher") << "hit strip csc from simhit" << endl;
0087 for (const auto& p : hit_strips) {
0088 edm::LogInfo("CSCRecHitMatcher") << p;
0089 }
0090 }
0091
0092
0093 const auto& rechits_in_det = rechits.get(p_id);
0094 for (auto d = rechits_in_det.first; d != rechits_in_det.second; ++d) {
0095 if (verboseCSCRecHit2D_)
0096 edm::LogInfo("CSCRecHitMatcher") << "rechit " << p_id << " " << *d;
0097
0098
0099 const bool wireMatch(std::find(hit_wg.begin(), hit_wg.end(), d->hitWire()) != hit_wg.end());
0100
0101
0102 bool stripMatch(false);
0103 for (size_t iS = 0; iS < d->nStrips(); ++iS) {
0104 if (hit_strips.find(d->channels(iS)) != hit_strips.end())
0105 stripMatch = true;
0106 }
0107
0108
0109 if (wireMatch and stripMatch) {
0110 if (verboseCSCRecHit2D_)
0111 edm::LogInfo("CSCRecHitMatcher") << "\t...was matched!";
0112 layer_to_cscRecHit2D_[id].push_back(*d);
0113 chamber_to_cscRecHit2D_[p_id.chamberId().rawId()].push_back(*d);
0114 }
0115 }
0116 }
0117 }
0118
0119 void CSCRecHitMatcher::matchCSCSegmentsToSimTrack(const CSCSegmentCollection& cscSegments) {
0120 if (verboseCSCSegment_)
0121 edm::LogInfo("CSCRecHitMatcher") << "Matching simtrack to segments";
0122
0123
0124 const auto& chamber_ids = chamberIdsCSCRecHit2D();
0125 if (verboseCSCSegment_)
0126 edm::LogInfo("CSCRecHitMatcher") << "Number of matched csc segments " << chamber_ids.size();
0127 for (const auto& id : chamber_ids) {
0128 CSCDetId p_id(id);
0129
0130
0131 const auto& csc_rechits(cscRecHit2DsInChamber(id));
0132 if (verboseCSCSegment_) {
0133 edm::LogInfo("CSCRecHitMatcher") << "hit csc rechits" << endl;
0134 for (const auto& p : csc_rechits) {
0135 edm::LogInfo("CSCRecHitMatcher") << p;
0136 }
0137 }
0138
0139
0140 const auto& segments_in_det = cscSegments.get(p_id);
0141 for (auto d = segments_in_det.first; d != segments_in_det.second; ++d) {
0142 if (verboseCSCSegment_)
0143 edm::LogInfo("CSCRecHitMatcher") << "segment " << p_id << " " << *d << endl;
0144
0145
0146 const auto& recHits(d->recHits());
0147
0148 int rechitsFound = 0;
0149 if (verboseCSCSegment_)
0150 edm::LogInfo("CSCRecHitMatcher") << recHits.size() << " csc rechits from segment " << endl;
0151 for (const auto& rh : recHits) {
0152 const CSCRecHit2D* cscrh(dynamic_cast<const CSCRecHit2D*>(rh));
0153 if (isCSCRecHit2DMatched(*cscrh))
0154 ++rechitsFound;
0155 }
0156 if (rechitsFound == 0)
0157 continue;
0158 if (verboseCSCSegment_) {
0159 edm::LogInfo("CSCRecHitMatcher") << "Found " << rechitsFound << " rechits out of "
0160 << cscRecHit2DsInChamber(id).size();
0161 edm::LogInfo("CSCRecHitMatcher") << "\t...was matched!";
0162 }
0163 chamber_to_cscSegment_[p_id.rawId()].push_back(*d);
0164 }
0165 }
0166 }
0167
0168 std::set<unsigned int> CSCRecHitMatcher::layerIdsCSCRecHit2D() const {
0169 std::set<unsigned int> result;
0170 for (const auto& p : layer_to_cscRecHit2D_)
0171 result.insert(p.first);
0172 return result;
0173 }
0174
0175 std::set<unsigned int> CSCRecHitMatcher::chamberIdsCSCRecHit2D() const {
0176 std::set<unsigned int> result;
0177 for (const auto& p : chamber_to_cscRecHit2D_)
0178 result.insert(p.first);
0179 return result;
0180 }
0181
0182 std::set<unsigned int> CSCRecHitMatcher::chamberIdsCSCSegment() const {
0183 std::set<unsigned int> result;
0184 for (const auto& p : chamber_to_cscSegment_)
0185 result.insert(p.first);
0186 return result;
0187 }
0188
0189 const CSCRecHit2DContainer& CSCRecHitMatcher::cscRecHit2DsInLayer(unsigned int detid) const {
0190 if (layer_to_cscRecHit2D_.find(detid) == layer_to_cscRecHit2D_.end())
0191 return no_cscRecHit2Ds_;
0192 return layer_to_cscRecHit2D_.at(detid);
0193 }
0194
0195 const CSCRecHit2DContainer& CSCRecHitMatcher::cscRecHit2DsInChamber(unsigned int detid) const {
0196 if (chamber_to_cscRecHit2D_.find(detid) == chamber_to_cscRecHit2D_.end())
0197 return no_cscRecHit2Ds_;
0198 return chamber_to_cscRecHit2D_.at(detid);
0199 }
0200
0201 const CSCSegmentContainer& CSCRecHitMatcher::cscSegmentsInChamber(unsigned int detid) const {
0202 if (chamber_to_cscSegment_.find(detid) == chamber_to_cscSegment_.end())
0203 return no_cscSegments_;
0204 return chamber_to_cscSegment_.at(detid);
0205 }
0206
0207 int CSCRecHitMatcher::nCSCRecHit2DsInLayer(unsigned int detid) const { return cscRecHit2DsInLayer(detid).size(); }
0208
0209 int CSCRecHitMatcher::nCSCRecHit2DsInChamber(unsigned int detid) const { return cscRecHit2DsInChamber(detid).size(); }
0210
0211 int CSCRecHitMatcher::nCSCSegmentsInChamber(unsigned int detid) const { return cscSegmentsInChamber(detid).size(); }
0212
0213 const CSCRecHit2DContainer CSCRecHitMatcher::cscRecHit2Ds() const {
0214 CSCRecHit2DContainer result;
0215 for (const auto& id : chamberIdsCSCRecHit2D()) {
0216 const auto& segmentsInChamber(cscRecHit2DsInChamber(id));
0217 result.insert(result.end(), segmentsInChamber.begin(), segmentsInChamber.end());
0218 }
0219 return result;
0220 }
0221
0222 const CSCSegmentContainer CSCRecHitMatcher::cscSegments() const {
0223 CSCSegmentContainer result;
0224 for (const auto& id : chamberIdsCSCSegment()) {
0225 const auto& segmentsInChamber(cscSegmentsInChamber(id));
0226 result.insert(result.end(), segmentsInChamber.begin(), segmentsInChamber.end());
0227 }
0228 return result;
0229 }
0230
0231 bool CSCRecHitMatcher::cscRecHit2DInContainer(const CSCRecHit2D& sg, const CSCRecHit2DContainer& c) const {
0232 bool isSame = false;
0233 for (const auto& segment : c)
0234 if (areCSCRecHit2DsSame(sg, segment))
0235 isSame = true;
0236 return isSame;
0237 }
0238
0239 bool CSCRecHitMatcher::cscSegmentInContainer(const CSCSegment& sg, const CSCSegmentContainer& c) const {
0240 bool isSame = false;
0241 for (const auto& segment : c)
0242 if (areCSCSegmentsSame(sg, segment))
0243 isSame = true;
0244 return isSame;
0245 }
0246
0247 bool CSCRecHitMatcher::isCSCRecHit2DMatched(const CSCRecHit2D& thisSg) const {
0248 return cscRecHit2DInContainer(thisSg, cscRecHit2Ds());
0249 }
0250
0251 bool CSCRecHitMatcher::isCSCSegmentMatched(const CSCSegment& thisSg) const {
0252 return cscSegmentInContainer(thisSg, cscSegments());
0253 }
0254
0255 int CSCRecHitMatcher::nCSCRecHit2Ds() const {
0256 int n = 0;
0257 const auto& ids = chamberIdsCSCRecHit2D();
0258 for (const auto& id : ids)
0259 n += cscRecHit2DsInChamber(id).size();
0260 return n;
0261 }
0262
0263 int CSCRecHitMatcher::nCSCSegments() const {
0264 int n = 0;
0265 const auto& ids = chamberIdsCSCSegment();
0266 for (const auto& id : ids)
0267 n += cscSegmentsInChamber(id).size();
0268 return n;
0269 }
0270
0271 bool CSCRecHitMatcher::areCSCRecHit2DsSame(const CSCRecHit2D& l, const CSCRecHit2D& r) const {
0272 return l.localPosition() == r.localPosition();
0273 }
0274
0275 bool CSCRecHitMatcher::areCSCSegmentsSame(const CSCSegment& l, const CSCSegment& r) const {
0276 return (l.localPosition() == r.localPosition() and l.localDirection() == r.localDirection());
0277 }
0278
0279 CSCSegment CSCRecHitMatcher::bestCSCSegment(unsigned int id) {
0280 CSCSegment emptySegment;
0281 double chi2overNdf = 99;
0282 int index = 0;
0283 int foundIndex = -99;
0284
0285 for (const auto& seg : chamber_to_cscSegment_[id]) {
0286 double newChi2overNdf(seg.chi2() / seg.degreesOfFreedom());
0287 if (newChi2overNdf < chi2overNdf) {
0288 chi2overNdf = newChi2overNdf;
0289 foundIndex = index;
0290 }
0291 ++index;
0292 }
0293 if (foundIndex != -99)
0294 return chamber_to_cscSegment_[id][foundIndex];
0295 return emptySegment;
0296 }
0297
0298 GlobalPoint CSCRecHitMatcher::globalPoint(const CSCSegment& c) const {
0299 return cscGeometry_->idToDet(c.cscDetId())->surface().toGlobal(c.localPosition());
0300 }