MuonCSCDigis/src/CSCDigiMatcher.cc

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