CSCRecHits/src/CSCRecHitMatcher.cc

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   // make a new digi matcher
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 /// do the matching
0042 void CSCRecHitMatcher::match(const SimTrack& t, const SimVertex& v) {
0043   // match digis first
0044   cscDigiMatcher_->match(t, v);
0045
0046   // get the rechit collection
0047   const CSCRecHit2DCollection& cscRecHit2Ds = *cscRecHit2DH_.product();
0048   const CSCSegmentCollection& cscSegments = *cscSegmentH_.product();
0049
0050   // now match the rechits
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   // fetch all layerIds with digis
0060   const auto& strip_ids = cscDigiMatcher_->detIdsStrip();
0061   const auto& wire_ids = cscDigiMatcher_->detIdsWire();
0062
0063   // merge the two collections
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     // print all the wires in the CSCChamber
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     // print all the strips in the CSCChamber
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     // get the rechits
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       // does the wire number match?
0099       const bool wireMatch(std::find(hit_wg.begin(), hit_wg.end(), d->hitWire()) != hit_wg.end());
0100
0101       // does the strip number match?
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       // this rechit was matched to a matching simhit
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   // fetch all chamberIds with 2D rechits
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     // print all CSCRecHit2D in the CSCChamber
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     // get the segments
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       //access the rechits
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 }