MuonCSCDigis/src/CSCDigiMatcher.cc

0001 #include <memory>
0002
0003 #include "Validation/MuonCSCDigis/interface/CSCDigiMatcher.h"
0004
0005 using namespace std;
0006
0007 CSCDigiMatcher::CSCDigiMatcher(const edm::ParameterSet& pset, edm::ConsumesCollector&& iC) {
0008   const auto& wireDigi = pset.getParameterSet("cscWireDigi");
0009   verboseWG_ = wireDigi.getParameter<int>("verbose");
0010   minBXWire_ = wireDigi.getParameter<int>("minBX");
0011   maxBXWire_ = wireDigi.getParameter<int>("maxBX");
0012   matchDeltaWG_ = wireDigi.getParameter<int>("matchDeltaWG");
0013
0014   const auto& comparatorDigi = pset.getParameterSet("cscComparatorDigi");
0015   verboseComparator_ = comparatorDigi.getParameter<int>("verbose");
0016   minBXComparator_ = comparatorDigi.getParameter<int>("minBX");
0017   maxBXComparator_ = comparatorDigi.getParameter<int>("maxBX");
0018   matchDeltaComparator_ = comparatorDigi.getParameter<int>("matchDeltaStrip");
0019
0020   const auto& stripDigi = pset.getParameterSet("cscStripDigi");
0021   verboseStrip_ = stripDigi.getParameter<int>("verbose");
0022   minBXStrip_ = stripDigi.getParameter<int>("minBX");
0023   maxBXStrip_ = stripDigi.getParameter<int>("maxBX");
0024   matchDeltaStrip_ = stripDigi.getParameter<int>("matchDeltaStrip");
0025
0026   // make a new simhits matcher
0027   muonSimHitMatcher_ = std::make_shared<CSCSimHitMatcher>(pset, std::move(iC));
0028
0029   comparatorDigiInput_ =
0030       iC.consumes<CSCComparatorDigiCollection>(comparatorDigi.getParameter<edm::InputTag>("inputTag"));
0031   stripDigiInput_ = iC.consumes<CSCStripDigiCollection>(stripDigi.getParameter<edm::InputTag>("inputTag"));
0032   wireDigiInput_ = iC.consumes<CSCWireDigiCollection>(wireDigi.getParameter<edm::InputTag>("inputTag"));
0033 }
0034
0035 void CSCDigiMatcher::init(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0036   muonSimHitMatcher_->init(iEvent, iSetup);
0037
0038   iEvent.getByToken(comparatorDigiInput_, comparatorDigisH_);
0039   iEvent.getByToken(stripDigiInput_, stripDigisH_);
0040   iEvent.getByToken(wireDigiInput_, wireDigisH_);
0041 }
0042
0043 /// do the matching
0044 void CSCDigiMatcher::match(const SimTrack& t, const SimVertex& v) {
0045   // match simhits first
0046   muonSimHitMatcher_->match(t, v);
0047
0048   // get the digi collections
0049   const CSCComparatorDigiCollection& comparators = *comparatorDigisH_.product();
0050   const CSCStripDigiCollection& strips = *stripDigisH_.product();
0051   const CSCWireDigiCollection& wires = *wireDigisH_.product();
0052
0053   clear();
0054
0055   // now match the digis
0056   matchComparatorsToSimTrack(comparators);
0057   matchStripsToSimTrack(strips);
0058   matchWiresToSimTrack(wires);
0059 }
0060
0061 void CSCDigiMatcher::matchComparatorsToSimTrack(const CSCComparatorDigiCollection& comparators) {
0062   const auto& det_ids = muonSimHitMatcher_->detIds(0);
0063   for (const auto& id : det_ids) {
0064     CSCDetId layer_id(id);
0065
0066     const auto& hit_comparators = muonSimHitMatcher_->hitStripsInDetId(id, matchDeltaStrip_);
0067     if (verboseComparator_) {
0068       cout << "hit_comparators_fat, CSCid " << layer_id << " ";
0069       copy(hit_comparators.begin(), hit_comparators.end(), ostream_iterator<int>(cout, " "));
0070       cout << endl;
0071     }
0072
0073     int ndigis = 0;
0074     const auto& comp_digis_in_det = comparators.get(layer_id);
0075     for (auto c = comp_digis_in_det.first; c != comp_digis_in_det.second; ++c) {
0076       if (verboseComparator_)
0077         edm::LogInfo("CSCDigiMatcher") << "sdigi " << layer_id << " (comparator, comparator, Tbin ) " << *c;
0078
0079       // check that the first BX for this digi wasn't too early or too late
0080       if (c->getTimeBin() < minBXComparator_ || c->getTimeBin() > maxBXComparator_)
0081         continue;
0082
0083       ndigis++;
0084
0085       int comparator = c->getStrip();  // comparators are counted from 1
0086       // check that it matches a comparator that was hit by SimHits from our track
0087       if (hit_comparators.find(comparator) == hit_comparators.end())
0088         continue;
0089
0090       if (verboseComparator_)
0091         edm::LogInfo("CSCDigiMatcher") << "Matched comparator " << *c;
0092       detid_to_comparators_[id].push_back(*c);
0093       chamber_to_comparators_[layer_id.chamberId().rawId()].push_back(*c);
0094     }
0095     detid_to_totalcomparators_[id] = ndigis;  //id to totalcomparators
0096   }
0097 }
0098
0099 void CSCDigiMatcher::matchStripsToSimTrack(const CSCStripDigiCollection& strips) {
0100   for (auto detUnitIt = strips.begin(); detUnitIt != strips.end(); ++detUnitIt) {
0101     const CSCDetId& id = (*detUnitIt).first;
0102     const auto& range = (*detUnitIt).second;
0103     for (auto digiIt = range.first; digiIt != range.second; ++digiIt) {
0104       if (id.station() == 1 and (id.ring() == 1 or id.ring() == 4))
0105         if (verboseStrip_)
0106           edm::LogInfo("CSCDigiMatcher") << "CSCid " << id << " Strip digi (strip, strip, Tbin ) " << (*digiIt);
0107     }
0108   }
0109
0110   const auto& det_ids = muonSimHitMatcher_->detIds(0);
0111   for (const auto& id : det_ids) {
0112     CSCDetId layer_id(id);
0113
0114     const auto& hit_strips = muonSimHitMatcher_->hitStripsInDetId(id, matchDeltaStrip_);
0115     if (verboseStrip_) {
0116       cout << "hit_strips_fat, CSCid " << layer_id << " ";
0117       copy(hit_strips.begin(), hit_strips.end(), ostream_iterator<int>(cout, " "));
0118       cout << endl;
0119     }
0120
0121     int ndigis = 0;
0122
0123     const auto& strip_digis_in_det = strips.get(layer_id);
0124     for (auto c = strip_digis_in_det.first; c != strip_digis_in_det.second; ++c) {
0125       //next is to remove the strips with pulse at noise level, with ACD info
0126       //the detail may be from CSC rechits algorithm
0127       if (verboseStrip_)
0128         edm::LogInfo("CSCDigiMatcher") << "sdigi " << layer_id << " (strip, ADC ) " << *c;
0129
0130       ndigis++;
0131
0132       int strip = c->getStrip();  // strips are counted from 1
0133       // check that it matches a strip that was hit by SimHits from our track
0134       if (hit_strips.find(strip) == hit_strips.end())
0135         continue;
0136
0137       if (verboseStrip_)
0138         edm::LogInfo("CSCDigiMatcher") << "Matched strip " << *c;
0139       detid_to_strips_[id].push_back(*c);
0140       chamber_to_strips_[layer_id.chamberId().rawId()].push_back(*c);
0141     }
0142     detid_to_totalstrips_[id] = ndigis;
0143   }
0144 }
0145
0146 void CSCDigiMatcher::matchWiresToSimTrack(const CSCWireDigiCollection& wires) {
0147   const auto& det_ids = muonSimHitMatcher_->detIds(0);
0148   for (const auto& id : det_ids) {
0149     CSCDetId layer_id(id);
0150
0151     const auto& hit_wires = muonSimHitMatcher_->hitWiregroupsInDetId(id, matchDeltaWG_);
0152     if (verboseWG_) {
0153       cout << "hit_wires ";
0154       copy(hit_wires.begin(), hit_wires.end(), ostream_iterator<int>(cout, " "));
0155       cout << endl;
0156     }
0157
0158     int ndigis = 0;
0159
0160     const auto& wire_digis_in_det = wires.get(layer_id);
0161     for (auto w = wire_digis_in_det.first; w != wire_digis_in_det.second; ++w) {
0162       if (verboseStrip_)
0163         edm::LogInfo("CSCDigiMatcher") << "wdigi " << layer_id << " (wire, Tbin ) " << *w;
0164
0165       // check that the first BX for this digi wasn't too early or too late
0166       if (w->getTimeBin() < minBXWire_ || w->getTimeBin() > maxBXWire_)
0167         continue;
0168
0169       ndigis++;
0170
0171       int wg = w->getWireGroup();  // wiregroups are counted from 1
0172       // check that it matches a strip that was hit by SimHits from our track
0173       if (hit_wires.find(wg) == hit_wires.end())
0174         continue;
0175
0176       if (verboseStrip_)
0177         edm::LogInfo("CSCDigiMatcher") << "Matched wire digi " << *w << endl;
0178       detid_to_wires_[id].push_back(*w);
0179       chamber_to_wires_[layer_id.chamberId().rawId()].push_back(*w);
0180     }
0181     detid_to_totalwires_[id] = ndigis;
0182   }
0183 }
0184
0185 std::set<unsigned int> CSCDigiMatcher::detIdsComparator(int csc_type) const {
0186   return selectDetIds(detid_to_comparators_, csc_type);
0187 }
0188
0189 std::set<unsigned int> CSCDigiMatcher::detIdsStrip(int csc_type) const {
0190   return selectDetIds(detid_to_strips_, csc_type);
0191 }
0192
0193 std::set<unsigned int> CSCDigiMatcher::detIdsWire(int csc_type) const {
0194   return selectDetIds(detid_to_wires_, csc_type);
0195 }
0196
0197 std::set<unsigned int> CSCDigiMatcher::chamberIdsComparator(int csc_type) const {
0198   return selectDetIds(chamber_to_comparators_, csc_type);
0199 }
0200
0201 std::set<unsigned int> CSCDigiMatcher::chamberIdsStrip(int csc_type) const {
0202   return selectDetIds(chamber_to_strips_, csc_type);
0203 }
0204
0205 std::set<unsigned int> CSCDigiMatcher::chamberIdsWire(int csc_type) const {
0206   return selectDetIds(chamber_to_wires_, csc_type);
0207 }
0208
0209 const CSCComparatorDigiContainer& CSCDigiMatcher::comparatorDigisInDetId(unsigned int detid) const {
0210   if (detid_to_comparators_.find(detid) == detid_to_comparators_.end())
0211     return no_comparators_;
0212   return detid_to_comparators_.at(detid);
0213 }
0214
0215 const CSCComparatorDigiContainer& CSCDigiMatcher::comparatorDigisInChamber(unsigned int detid) const {
0216   if (chamber_to_comparators_.find(detid) == chamber_to_comparators_.end())
0217     return no_comparators_;
0218   return chamber_to_comparators_.at(detid);
0219 }
0220
0221 const CSCStripDigiContainer& CSCDigiMatcher::stripDigisInDetId(unsigned int detid) const {
0222   if (detid_to_strips_.find(detid) == detid_to_strips_.end())
0223     return no_strips_;
0224   return detid_to_strips_.at(detid);
0225 }
0226
0227 const CSCStripDigiContainer& CSCDigiMatcher::stripDigisInChamber(unsigned int detid) const {
0228   if (chamber_to_strips_.find(detid) == chamber_to_strips_.end())
0229     return no_strips_;
0230   return chamber_to_strips_.at(detid);
0231 }
0232
0233 const CSCWireDigiContainer& CSCDigiMatcher::wireDigisInDetId(unsigned int detid) const {
0234   if (detid_to_wires_.find(detid) == detid_to_wires_.end())
0235     return no_wires_;
0236   return detid_to_wires_.at(detid);
0237 }
0238
0239 const CSCWireDigiContainer& CSCDigiMatcher::wireDigisInChamber(unsigned int detid) const {
0240   if (chamber_to_wires_.find(detid) == chamber_to_wires_.end())
0241     return no_wires_;
0242   return chamber_to_wires_.at(detid);
0243 }
0244
0245 int CSCDigiMatcher::nLayersWithComparatorInChamber(unsigned int detid) const {
0246   int nLayers = 0;
0247   CSCDetId chamberId(detid);
0248   for (int i = 1; i <= 6; ++i) {
0249     CSCDetId layerId(chamberId.endcap(), chamberId.station(), chamberId.ring(), chamberId.chamber(), i);
0250     if (!comparatorDigisInDetId(layerId.rawId()).empty()) {
0251       nLayers++;
0252     }
0253   }
0254   return nLayers;
0255 }
0256
0257 int CSCDigiMatcher::nLayersWithStripInChamber(unsigned int detid) const {
0258   int nLayers = 0;
0259   CSCDetId chamberId(detid);
0260   for (int i = 1; i <= 6; ++i) {
0261     CSCDetId layerId(chamberId.endcap(), chamberId.station(), chamberId.ring(), chamberId.chamber(), i);
0262     if (!stripDigisInDetId(layerId.rawId()).empty()) {
0263       nLayers++;
0264     }
0265   }
0266   return nLayers;
0267 }
0268
0269 int CSCDigiMatcher::nLayersWithWireInChamber(unsigned int detid) const {
0270   int nLayers = 0;
0271   CSCDetId chamberId(detid);
0272   for (int i = 1; i <= 6; ++i) {
0273     CSCDetId layerId(chamberId.endcap(), chamberId.station(), chamberId.ring(), chamberId.chamber(), i);
0274     if (!wireDigisInDetId(layerId.rawId()).empty()) {
0275       nLayers++;
0276     }
0277   }
0278   return nLayers;
0279 }
0280
0281 int CSCDigiMatcher::nCoincidenceComparatorChambers(int min_n_layers) const {
0282   int result = 0;
0283   const auto& chamber_ids = chamberIdsComparator();
0284   for (const auto& id : chamber_ids) {
0285     if (nLayersWithComparatorInChamber(id) >= min_n_layers)
0286       result += 1;
0287   }
0288   return result;
0289 }
0290
0291 int CSCDigiMatcher::nCoincidenceStripChambers(int min_n_layers) const {
0292   int result = 0;
0293   const auto& chamber_ids = chamberIdsStrip();
0294   for (const auto& id : chamber_ids) {
0295     if (nLayersWithStripInChamber(id) >= min_n_layers)
0296       result += 1;
0297   }
0298   return result;
0299 }
0300
0301 int CSCDigiMatcher::nCoincidenceWireChambers(int min_n_layers) const {
0302   int result = 0;
0303   const auto& chamber_ids = chamberIdsWire();
0304   for (const auto& id : chamber_ids) {
0305     if (nLayersWithWireInChamber(id) >= min_n_layers)
0306       result += 1;
0307   }
0308   return result;
0309 }
0310
0311 std::set<int> CSCDigiMatcher::comparatorsInDetId(unsigned int detid) const {
0312   set<int> result;
0313   const auto& digis = comparatorDigisInDetId(detid);
0314   for (const auto& d : digis) {
0315     result.insert(d.getHalfStrip());
0316   }
0317   return result;
0318 }
0319
0320 std::set<int> CSCDigiMatcher::stripsInDetId(unsigned int detid) const {
0321   set<int> result;
0322   const auto& digis = stripDigisInDetId(detid);
0323   for (const auto& d : digis) {
0324     result.insert(d.getStrip());
0325   }
0326   return result;
0327 }
0328
0329 std::set<int> CSCDigiMatcher::wiregroupsInDetId(unsigned int detid) const {
0330   set<int> result;
0331   const auto& digis = wireDigisInDetId(detid);
0332   for (const auto& d : digis) {
0333     result.insert(d.getWireGroup());
0334   }
0335   return result;
0336 }
0337
0338 std::set<int> CSCDigiMatcher::comparatorsInChamber(unsigned int detid, int max_gap_to_fill) const {
0339   set<int> result;
0340   const auto& digis = comparatorDigisInChamber(detid);
0341   for (const auto& d : digis) {
0342     result.insert(d.getStrip());
0343   }
0344   if (max_gap_to_fill > 0) {
0345     int prev = -111;
0346     for (const auto& s : result) {
0347       if (s - prev > 1 && s - prev - 1 <= max_gap_to_fill) {
0348         for (int fill_s = prev + 1; fill_s < s; ++fill_s)
0349           result.insert(fill_s);
0350       }
0351       prev = s;
0352     }
0353   }
0354
0355   return result;
0356 }
0357
0358 std::set<int> CSCDigiMatcher::stripsInChamber(unsigned int detid, int max_gap_to_fill) const {
0359   set<int> result;
0360   const auto& digis = stripDigisInChamber(detid);
0361   for (const auto& d : digis) {
0362     result.insert(d.getStrip());
0363   }
0364   if (max_gap_to_fill > 0) {
0365     int prev = -111;
0366     for (const auto& s : result) {
0367       if (s - prev > 1 && s - prev - 1 <= max_gap_to_fill) {
0368         for (int fill_s = prev + 1; fill_s < s; ++fill_s)
0369           result.insert(fill_s);
0370       }
0371       prev = s;
0372     }
0373   }
0374
0375   return result;
0376 }
0377
0378 std::set<int> CSCDigiMatcher::wiregroupsInChamber(unsigned int detid, int max_gap_to_fill) const {
0379   set<int> result;
0380   const auto& digis = wireDigisInChamber(detid);
0381   for (const auto& d : digis) {
0382     result.insert(d.getWireGroup());
0383   }
0384   if (max_gap_to_fill > 0) {
0385     int prev = -111;
0386     for (const auto& w : result) {
0387       if (w - prev > 1 && w - prev - 1 <= max_gap_to_fill) {
0388         for (int fill_w = prev + 1; fill_w < w; ++fill_w)
0389           result.insert(fill_w);
0390       }
0391       prev = w;
0392     }
0393   }
0394   return result;
0395 }
0396
0397 int CSCDigiMatcher::totalComparators(unsigned int detid) const {
0398   if (detid_to_totalcomparators_.find(detid) == detid_to_totalcomparators_.end())
0399     return 0;
0400   return detid_to_totalcomparators_.at(detid);
0401 }
0402
0403 int CSCDigiMatcher::totalStrips(unsigned int detid) const {
0404   if (detid_to_totalstrips_.find(detid) == detid_to_totalstrips_.end())
0405     return 0;
0406   return detid_to_totalstrips_.at(detid);
0407 }
0408
0409 int CSCDigiMatcher::totalWires(unsigned int detid) const {
0410   if (detid_to_totalwires_.find(detid) == detid_to_totalwires_.end())
0411     return 0;
0412   return detid_to_totalwires_.at(detid);
0413 }
0414
0415 void CSCDigiMatcher::clear() {
0416   detid_to_comparators_.clear();
0417   chamber_to_comparators_.clear();
0418
0419   detid_to_strips_.clear();
0420   chamber_to_strips_.clear();
0421
0422   detid_to_wires_.clear();
0423   chamber_to_wires_.clear();
0424 }