MuonCSCDigis/src/CSCStubMatcher.cc

0001 #include "Validation/MuonCSCDigis/interface/CSCStubMatcher.h"
0002 #include "DataFormats/CSCDigi/interface/CSCConstants.h"
0003 #include <algorithm>
0004
0005 using namespace std;
0006
0007 CSCStubMatcher::CSCStubMatcher(const edm::ParameterSet& pSet, edm::ConsumesCollector&& iC) {
0008   useGEMs_ = pSet.getParameter<bool>("useGEMs");
0009
0010   const auto& cscCLCT = pSet.getParameter<edm::ParameterSet>("cscCLCT");
0011   minBXCLCT_ = cscCLCT.getParameter<int>("minBX");
0012   maxBXCLCT_ = cscCLCT.getParameter<int>("maxBX");
0013   verboseCLCT_ = cscCLCT.getParameter<int>("verbose");
0014   minNHitsChamberCLCT_ = cscCLCT.getParameter<int>("minNHitsChamber");
0015
0016   const auto& cscALCT = pSet.getParameter<edm::ParameterSet>("cscALCT");
0017   minBXALCT_ = cscALCT.getParameter<int>("minBX");
0018   maxBXALCT_ = cscALCT.getParameter<int>("maxBX");
0019   verboseALCT_ = cscALCT.getParameter<int>("verbose");
0020   minNHitsChamberALCT_ = cscALCT.getParameter<int>("minNHitsChamber");
0021
0022   const auto& cscLCT = pSet.getParameter<edm::ParameterSet>("cscLCT");
0023   minBXLCT_ = cscLCT.getParameter<int>("minBX");
0024   maxBXLCT_ = cscLCT.getParameter<int>("maxBX");
0025   verboseLCT_ = cscLCT.getParameter<int>("verbose");
0026   minNHitsChamberLCT_ = cscLCT.getParameter<int>("minNHitsChamber");
0027   addGhostLCTs_ = cscLCT.getParameter<bool>("addGhosts");
0028
0029   const auto& cscMPLCT = pSet.getParameter<edm::ParameterSet>("cscMPLCT");
0030   minBXMPLCT_ = cscMPLCT.getParameter<int>("minBX");
0031   maxBXMPLCT_ = cscMPLCT.getParameter<int>("maxBX");
0032   verboseMPLCT_ = cscMPLCT.getParameter<int>("verbose");
0033   minNHitsChamberMPLCT_ = cscMPLCT.getParameter<int>("minNHitsChamber");
0034
0035   if (useGEMs_)
0036     gemDigiMatcher_.reset(new GEMDigiMatcher(pSet, std::move(iC)));
0037   cscDigiMatcher_.reset(new CSCDigiMatcher(pSet, std::move(iC)));
0038
0039   clctInputTag_ = cscCLCT.getParameter<edm::InputTag>("inputTag");
0040   alctInputTag_ = cscALCT.getParameter<edm::InputTag>("inputTag");
0041   lctInputTag_ = cscLCT.getParameter<edm::InputTag>("inputTag");
0042   mplctInputTag_ = cscMPLCT.getParameter<edm::InputTag>("inputTag");
0043
0044   clctToken_ = iC.consumes<CSCCLCTDigiCollection>(clctInputTag_);
0045   alctToken_ = iC.consumes<CSCALCTDigiCollection>(alctInputTag_);
0046   lctToken_ = iC.consumes<CSCCorrelatedLCTDigiCollection>(lctInputTag_);
0047   mplctToken_ = iC.consumes<CSCCorrelatedLCTDigiCollection>(mplctInputTag_);
0048
0049   geomToken_ = iC.esConsumes<CSCGeometry, MuonGeometryRecord>();
0050 }
0051
0052 void CSCStubMatcher::init(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0053   if (useGEMs_)
0054     gemDigiMatcher_->init(iEvent, iSetup);
0055   cscDigiMatcher_->init(iEvent, iSetup);
0056
0057   iEvent.getByToken(clctToken_, clctsH_);
0058   iEvent.getByToken(alctToken_, alctsH_);
0059   iEvent.getByToken(lctToken_, lctsH_);
0060   iEvent.getByToken(mplctToken_, mplctsH_);
0061
0062   cscGeometry_ = &iSetup.getData(geomToken_);
0063 }
0064
0065 // do the matching
0066 void CSCStubMatcher::match(const SimTrack& t, const SimVertex& v) {
0067   // match simhits first
0068   if (useGEMs_)
0069     gemDigiMatcher_->match(t, v);
0070   cscDigiMatcher_->match(t, v);
0071
0072   const CSCCLCTDigiCollection& clcts = *clctsH_.product();
0073   const CSCALCTDigiCollection& alcts = *alctsH_.product();
0074   const CSCCorrelatedLCTDigiCollection& lcts = *lctsH_.product();
0075   const CSCCorrelatedLCTDigiCollection& mplcts = *mplctsH_.product();
0076
0077   // clear collections
0078   clear();
0079
0080   if (!alctsH_.isValid()) {
0081     edm::LogError("CSCStubMatcher") << "Cannot get ALCTs with label " << alctInputTag_.encode();
0082   } else {
0083     matchALCTsToSimTrack(alcts);
0084   }
0085
0086   if (!clctsH_.isValid()) {
0087     edm::LogError("CSCStubMatcher") << "Cannot get CLCTs with label " << clctInputTag_.encode();
0088   } else {
0089     matchCLCTsToSimTrack(clcts);
0090   }
0091
0092   if (!lctsH_.isValid()) {
0093     edm::LogError("CSCStubMatcher") << "Cannot get LCTs with label " << lctInputTag_.encode();
0094   } else {
0095     matchLCTsToSimTrack(lcts);
0096   }
0097
0098   if (!mplctsH_.isValid()) {
0099     edm::LogError("CSCStubMatcher") << "Cannot get MPLCTs with label " << mplctInputTag_.encode();
0100   } else {
0101     matchMPLCTsToSimTrack(mplcts);
0102   }
0103 }
0104
0105 void CSCStubMatcher::matchCLCTsToSimTrack(const CSCCLCTDigiCollection& clcts) {
0106   const auto& cathode_ids = cscDigiMatcher_->chamberIdsStrip(0);
0107
0108   for (const auto& id : cathode_ids) {
0109     CSCDetId ch_id(id);
0110     if (verboseCLCT_) {
0111       edm::LogInfo("CSCStubMatcher") << "To check CSC chamber " << ch_id;
0112     }
0113
0114     int ring = ch_id.ring();
0115
0116     // do not consider CSCs with too few hits
0117     if (cscDigiMatcher_->nLayersWithStripInChamber(ch_id) < minNHitsChamberCLCT_)
0118       continue;
0119
0120     // get the comparator digis in this chamber
0121     std::vector<CSCComparatorDigiContainer> comps;
0122     for (int ilayer = CSCDetId::minLayerId(); ilayer <= CSCDetId::maxLayerId(); ilayer++) {
0123       CSCDetId layerid(ch_id.endcap(), ch_id.station(), ring, ch_id.chamber(), ilayer);
0124       comps.push_back(cscDigiMatcher_->comparatorDigisInDetId(layerid));
0125     }
0126
0127     // print out the digis
0128     if (verboseCLCT_) {
0129       edm::LogInfo("CSCStubMatcher") << "clct: comparators " << ch_id;
0130       int layer = 0;
0131       for (const auto& p : comps) {
0132         layer++;
0133         for (const auto& q : p) {
0134           edm::LogInfo("CSCStubMatcher") << "L" << layer << " " << q << " " << q.getHalfStrip() << " ";
0135         }
0136       }
0137     }
0138
0139     //use ME1b id to get CLCTs
0140     const bool isME1a(ch_id.station() == 1 and ch_id.ring() == 4);
0141     if (isME1a)
0142       ring = 1;
0143     CSCDetId ch_id2(ch_id.endcap(), ch_id.station(), ring, ch_id.chamber(), 0);
0144     auto id2 = ch_id2.rawId();  // CLCTs should be sorted into the det of the CLCTs.
0145
0146     const auto& clcts_in_det = clcts.get(ch_id2);
0147
0148     for (auto c = clcts_in_det.first; c != clcts_in_det.second; ++c) {
0149       if (verboseCLCT_)
0150         edm::LogInfo("CSCStubMatcher") << "clct " << ch_id2 << " " << *c;
0151
0152       if (!c->isValid())
0153         continue;
0154
0155       // check that the BX for this stub wasn't too early or too late
0156       if (c->getBX() < minBXCLCT_ || c->getBX() > maxBXCLCT_)
0157         continue;
0158
0159       // store all CLCTs in this chamber
0160       chamber_to_clcts_all_[id2].push_back(*c);
0161
0162       // check that at least 3 comparator digis were matched!
0163       int nMatches = 0;
0164       int layer = 0;
0165       for (const auto& p : comps) {
0166         layer++;
0167         for (const auto& q : p) {
0168           if (verboseCLCT_)
0169             edm::LogInfo("CSCStubMatcher") << "L" << layer << " " << q << " " << q.getHalfStrip() << " " << std::endl;
0170           for (const auto& clctComp : (*c).getHits()[layer - 1]) {
0171             if (clctComp == 65535)
0172               continue;
0173             if (verboseCLCT_) {
0174               edm::LogInfo("CSCStubMatcher") << "\t" << clctComp << " ";
0175             }
0176             if (q.getHalfStrip() == clctComp or (isME1a and q.getHalfStrip() + 128 == clctComp)) {
0177               nMatches++;
0178               if (verboseCLCT_) {
0179                 edm::LogInfo("CSCStubMatcher") << "\t\tnMatches " << nMatches << std::endl;
0180               }
0181             }
0182           }
0183         }
0184       }
0185
0186       // require at least 3 good matches
0187       if (nMatches < 3)
0188         continue;
0189
0190       if (verboseCLCT_)
0191         edm::LogInfo("CSCStubMatcher") << "clctGOOD";
0192
0193       // store matching CLCTs in this chamber
0194       if (std::find(chamber_to_clcts_[id2].begin(), chamber_to_clcts_[id2].end(), *c) == chamber_to_clcts_[id2].end()) {
0195         chamber_to_clcts_[id2].push_back(*c);
0196       }
0197     }
0198     if (chamber_to_clcts_[id2].size() > 2) {
0199       edm::LogInfo("CSCStubMatcher") << "WARNING!!! too many CLCTs " << chamber_to_clcts_[id2].size() << " in "
0200                                      << ch_id2;
0201       for (auto& c : chamber_to_clcts_[id2])
0202         edm::LogInfo("CSCStubMatcher") << "  " << c;
0203     }
0204   }
0205 }
0206
0207 void CSCStubMatcher::matchALCTsToSimTrack(const CSCALCTDigiCollection& alcts) {
0208   const auto& anode_ids = cscDigiMatcher_->chamberIdsWire(0);
0209   int n_minLayers = 0;
0210   for (const auto& id : anode_ids) {
0211     if (cscDigiMatcher_->nLayersWithWireInChamber(id) >= minNHitsChamberALCT_)
0212       ++n_minLayers;
0213     CSCDetId ch_id(id);
0214
0215     // fill 1 WG wide gaps
0216     const auto& digi_wgs = cscDigiMatcher_->wiregroupsInChamber(id, 1);
0217     if (verboseALCT_) {
0218       cout << "alct: digi_wgs " << ch_id << " ";
0219       copy(digi_wgs.begin(), digi_wgs.end(), ostream_iterator<int>(cout, " "));
0220       cout << endl;
0221     }
0222
0223     int ring = ch_id.ring();
0224     if (ring == 4)
0225       ring = 1;  //use ME1b id to get ALCTs
0226     CSCDetId ch_id2(ch_id.endcap(), ch_id.station(), ring, ch_id.chamber(), 0);
0227     auto id2 = ch_id2.rawId();  // ALCTs should be sorted into the det of the ALCTs.
0228
0229     const auto& alcts_in_det = alcts.get(ch_id2);
0230     for (auto a = alcts_in_det.first; a != alcts_in_det.second; ++a) {
0231       if (!a->isValid())
0232         continue;
0233
0234       if (verboseALCT_)
0235         edm::LogInfo("CSCStubMatcher") << "alct " << ch_id << " " << *a;
0236
0237       // check that the BX for stub wasn't too early or too late
0238       if (a->getBX() < minBXALCT_ || a->getBX() > maxBXALCT_)
0239         continue;
0240
0241       int wg = a->getKeyWG() + 1;  // as ALCT wiregroups numbers start from 0
0242
0243       // store all ALCTs in this chamber
0244       chamber_to_alcts_all_[id2].push_back(*a);
0245
0246       // match by wiregroup with the digis
0247       if (digi_wgs.find(wg) == digi_wgs.end()) {
0248         continue;
0249       }
0250       if (verboseALCT_)
0251         edm::LogInfo("CSCStubMatcher") << "alctGOOD";
0252
0253       // store matching ALCTs in this chamber
0254       if (std::find(chamber_to_alcts_[id2].begin(), chamber_to_alcts_[id2].end(), *a) == chamber_to_alcts_[id2].end()) {
0255         chamber_to_alcts_[id2].push_back(*a);
0256       }
0257     }
0258     if (chamber_to_alcts_[id2].size() > 2) {
0259       edm::LogInfo("CSCStubMatcher") << "WARNING!!! too many ALCTs " << chamber_to_alcts_[id2].size() << " in "
0260                                      << ch_id;
0261       for (auto& a : chamber_to_alcts_[id2])
0262         edm::LogInfo("CSCStubMatcher") << "  " << a;
0263     }
0264   }
0265 }
0266
0267 void CSCStubMatcher::matchLCTsToSimTrack(const CSCCorrelatedLCTDigiCollection& lcts) {
0268   // only look for stubs in chambers that already have CLCT and ALCT
0269   const auto& cathode_ids = chamberIdsAllCLCT(0);
0270   const auto& anode_ids = chamberIdsAllALCT(0);
0271
0272   std::set<int> cathode_and_anode_ids;
0273   std::set_union(cathode_ids.begin(),
0274                  cathode_ids.end(),
0275                  anode_ids.begin(),
0276                  anode_ids.end(),
0277                  std::inserter(cathode_and_anode_ids, cathode_and_anode_ids.end()));
0278
0279   for (const auto& id : cathode_and_anode_ids) {
0280     int iLct = -1;
0281
0282     CSCDetId ch_id(id);
0283
0284     //use ME1b id to get LCTs
0285     int ring = ch_id.ring();
0286     if (ring == 4)
0287       ring = 1;
0288     CSCDetId ch_id2(ch_id.endcap(), ch_id.station(), ring, ch_id.chamber(), 0);
0289     auto id2 = ch_id2.rawId();  // LCTs should be sorted into the det of the LCTs.
0290
0291     const auto& lcts_in_det = lcts.get(ch_id2);
0292
0293     std::map<int, CSCCorrelatedLCTDigiContainer> bx_to_lcts;
0294
0295     // collect all valid LCTs in a handy container
0296     CSCCorrelatedLCTDigiContainer lcts_tmp;
0297     for (auto lct = lcts_in_det.first; lct != lcts_in_det.second; ++lct) {
0298       if (!lct->isValid())
0299         continue;
0300       lcts_tmp.push_back(*lct);
0301       int bx = lct->getBX();
0302       bx_to_lcts[bx].push_back(*lct);
0303
0304       // Add ghost LCTs when there are two in bx
0305       // and the two don't share half-strip or wiregroup
0306       if (bx_to_lcts[bx].size() == 2 and addGhostLCTs_) {
0307         auto lct11 = bx_to_lcts[bx][0];
0308         auto lct22 = bx_to_lcts[bx][1];
0309         addGhostLCTs(lct11, lct22, lcts_tmp);
0310       }
0311     }
0312
0313     for (const auto& lct : lcts_tmp) {
0314       iLct++;
0315
0316       bool lct_clct_match(false);
0317       bool lct_alct_match(false);
0318       bool lct_gem1_match(false);
0319       bool lct_gem2_match(false);
0320
0321       if (verboseLCT_) {
0322         edm::LogInfo("CSCStubMatcher") << ch_id << " " << ch_id2;
0323         edm::LogInfo("CSCStubMatcher") << lct;
0324         edm::LogInfo("CSCStubMatcher") << "getCLCT " << lct.getCLCT() << "\ngetALCT " << lct.getALCT() << "\ngetGEM1 "
0325                                        << lct.getGEM1() << "\ngetGEM2 " << lct.getGEM2();
0326       }
0327       // Check if matched to an CLCT
0328       for (const auto& p : clctsInChamber(id)) {
0329         if (p == lct.getCLCT()) {
0330           lct_clct_match = true;
0331           if (verboseLCT_)
0332             edm::LogInfo("CSCStubMatcher") << "\t...lct_clct_match";
0333           break;
0334         }
0335       }
0336
0337       // Check if matched to an ALCT
0338       for (const auto& p : alctsInChamber(id)) {
0339         if (p == lct.getALCT()) {
0340           lct_alct_match = true;
0341           if (verboseLCT_)
0342             edm::LogInfo("CSCStubMatcher") << "\t...lct_alct_match";
0343           break;
0344         }
0345       }
0346
0347       if (useGEMs_) {
0348         // fixME here: double check the timing of GEMPad
0349         if (ch_id.ring() == 1 and (ch_id.station() == 1 or ch_id.station() == 2)) {
0350           // Check if matched to an GEM pad L1
0351           const GEMDetId gemDetIdL1(ch_id.zendcap(), 1, ch_id.station(), 1, ch_id.chamber(), 0);
0352           for (const auto& p : gemDigiMatcher_->padsInChamber(gemDetIdL1.rawId())) {
0353             if (p == lct.getGEM1()) {
0354               lct_gem1_match = true;
0355               if (verboseLCT_)
0356                 edm::LogInfo("CSCStubMatcher") << "\t...lct_gem1_match";
0357               break;
0358             }
0359           }
0360
0361           // Check if matched to an GEM pad L2
0362           const GEMDetId gemDetIdL2(ch_id.zendcap(), 1, ch_id.station(), 2, ch_id.chamber(), 0);
0363           for (const auto& p : gemDigiMatcher_->padsInChamber(gemDetIdL2.rawId())) {
0364             if (p == lct.getGEM2()) {
0365               lct_gem2_match = true;
0366               if (verboseLCT_)
0367                 edm::LogInfo("CSCStubMatcher") << "\t...lct_gem2_match";
0368               break;
0369             }
0370           }
0371         }
0372       }
0373
0374       const bool alct_clct = lct_clct_match and lct_alct_match;
0375       const bool alct_gem = lct_alct_match and lct_gem1_match and lct_gem2_match;
0376       const bool clct_gem = lct_clct_match and lct_gem1_match and lct_gem2_match;
0377
0378       bool lct_tight_matched = alct_clct or alct_gem or clct_gem;
0379       bool lct_loose_matched = lct_clct_match or lct_alct_match;
0380       bool lct_matched = lct_loose_matched or lct_tight_matched;
0381
0382       if (lct_matched) {
0383         if (verboseLCT_)
0384           edm::LogInfo("CSCStubMatcher") << "...was matched";
0385         if (std::find(chamber_to_lcts_[id2].begin(), chamber_to_lcts_[id2].end(), lct) == chamber_to_lcts_[id2].end()) {
0386           chamber_to_lcts_[id2].emplace_back(lct);
0387         }
0388       }
0389     }  // lct loop over
0390   }
0391 }
0392
0393 void CSCStubMatcher::matchMPLCTsToSimTrack(const CSCCorrelatedLCTDigiCollection& mplcts) {
0394   // match simtrack to MPC LCT by looking only in chambers
0395   // that already have LCTs matched to this simtrack
0396   const auto& lcts_ids = chamberIdsLCT(0);
0397
0398   // loop on the detids
0399   for (const auto& id : lcts_ids) {
0400     const auto& mplcts_in_det = mplcts.get(id);
0401
0402     // loop on the MPC LCTs in this detid
0403     for (auto lct = mplcts_in_det.first; lct != mplcts_in_det.second; ++lct) {
0404       if (!lct->isValid())
0405         continue;
0406
0407       chamber_to_mplcts_all_[id].emplace_back(*lct);
0408
0409       // check if this stub corresponds with a previously matched stub
0410       for (const auto& sim_stub : lctsInChamber(id)) {
0411         if (sim_stub == *lct) {
0412           if (std::find(chamber_to_mplcts_[id].begin(), chamber_to_mplcts_[id].end(), *lct) ==
0413               chamber_to_mplcts_[id].end()) {
0414             chamber_to_mplcts_[id].emplace_back(*lct);
0415           }
0416         }
0417       }
0418     }
0419   }
0420 }
0421
0422 std::set<unsigned int> CSCStubMatcher::chamberIdsAllCLCT(int csc_type) const {
0423   return selectDetIds(chamber_to_clcts_all_, csc_type);
0424 }
0425
0426 std::set<unsigned int> CSCStubMatcher::chamberIdsAllALCT(int csc_type) const {
0427   return selectDetIds(chamber_to_alcts_all_, csc_type);
0428 }
0429
0430 std::set<unsigned int> CSCStubMatcher::chamberIdsAllLCT(int csc_type) const {
0431   return selectDetIds(chamber_to_lcts_all_, csc_type);
0432 }
0433
0434 std::set<unsigned int> CSCStubMatcher::chamberIdsAllMPLCT(int csc_type) const {
0435   return selectDetIds(chamber_to_mplcts_all_, csc_type);
0436 }
0437
0438 std::set<unsigned int> CSCStubMatcher::chamberIdsCLCT(int csc_type) const {
0439   return selectDetIds(chamber_to_clcts_, csc_type);
0440 }
0441
0442 std::set<unsigned int> CSCStubMatcher::chamberIdsALCT(int csc_type) const {
0443   return selectDetIds(chamber_to_alcts_, csc_type);
0444 }
0445
0446 std::set<unsigned int> CSCStubMatcher::chamberIdsLCT(int csc_type) const {
0447   return selectDetIds(chamber_to_lcts_, csc_type);
0448 }
0449
0450 std::set<unsigned int> CSCStubMatcher::chamberIdsMPLCT(int csc_type) const {
0451   return selectDetIds(chamber_to_mplcts_, csc_type);
0452 }
0453
0454 const CSCCLCTDigiContainer& CSCStubMatcher::allCLCTsInChamber(unsigned int detid) const {
0455   if (chamber_to_clcts_all_.find(detid) == chamber_to_clcts_all_.end())
0456     return no_clcts_;
0457   return chamber_to_clcts_all_.at(detid);
0458 }
0459
0460 const CSCALCTDigiContainer& CSCStubMatcher::allALCTsInChamber(unsigned int detid) const {
0461   if (chamber_to_alcts_all_.find(detid) == chamber_to_alcts_all_.end())
0462     return no_alcts_;
0463   return chamber_to_alcts_all_.at(detid);
0464 }
0465
0466 const CSCCorrelatedLCTDigiContainer& CSCStubMatcher::allLCTsInChamber(unsigned int detid) const {
0467   if (chamber_to_lcts_all_.find(detid) == chamber_to_lcts_all_.end())
0468     return no_lcts_;
0469   return chamber_to_lcts_all_.at(detid);
0470 }
0471
0472 const CSCCorrelatedLCTDigiContainer& CSCStubMatcher::allMPLCTsInChamber(unsigned int detid) const {
0473   if (chamber_to_mplcts_all_.find(detid) == chamber_to_mplcts_all_.end())
0474     return no_mplcts_;
0475   return chamber_to_mplcts_all_.at(detid);
0476 }
0477
0478 const CSCCLCTDigiContainer& CSCStubMatcher::clctsInChamber(unsigned int detid) const {
0479   if (chamber_to_clcts_.find(detid) == chamber_to_clcts_.end())
0480     return no_clcts_;
0481   return chamber_to_clcts_.at(detid);
0482 }
0483
0484 const CSCALCTDigiContainer& CSCStubMatcher::alctsInChamber(unsigned int detid) const {
0485   if (chamber_to_alcts_.find(detid) == chamber_to_alcts_.end())
0486     return no_alcts_;
0487   return chamber_to_alcts_.at(detid);
0488 }
0489
0490 const CSCCorrelatedLCTDigiContainer& CSCStubMatcher::lctsInChamber(unsigned int detid) const {
0491   if (chamber_to_lcts_.find(detid) == chamber_to_lcts_.end())
0492     return no_lcts_;
0493   return chamber_to_lcts_.at(detid);
0494 }
0495
0496 const CSCCorrelatedLCTDigiContainer& CSCStubMatcher::mplctsInChamber(unsigned int detid) const {
0497   if (chamber_to_mplcts_.find(detid) == chamber_to_mplcts_.end())
0498     return no_mplcts_;
0499   return chamber_to_mplcts_.at(detid);
0500 }
0501
0502 CSCCLCTDigi CSCStubMatcher::bestClctInChamber(unsigned int detid) const {
0503   //sort stubs based on quality
0504   const auto& input(clctsInChamber(detid));
0505   int bestQ = 0;
0506   int index = -1;
0507   for (unsigned int i = 0; i < input.size(); ++i) {
0508     int quality = input[i].getQuality();
0509     if (quality > bestQ) {
0510       bestQ = quality;
0511       index = i;
0512     }
0513   }
0514   if (index != -1)
0515     return input[index];
0516   return CSCCLCTDigi();
0517 }
0518
0519 CSCALCTDigi CSCStubMatcher::bestAlctInChamber(unsigned int detid) const {
0520   //sort stubs based on quality
0521   const auto& input(alctsInChamber(detid));
0522   int bestQ = 0;
0523   int index = -1;
0524   for (unsigned int i = 0; i < input.size(); ++i) {
0525     int quality = input[i].getQuality();
0526     if (quality > bestQ) {
0527       bestQ = quality;
0528       index = i;
0529     }
0530   }
0531   if (index != -1)
0532     return input[index];
0533   return CSCALCTDigi();
0534 }
0535
0536 CSCCorrelatedLCTDigi CSCStubMatcher::bestLctInChamber(unsigned int detid) const {
0537   //sort stubs based on quality
0538   const auto& input(lctsInChamber(detid));
0539   int bestQ = 0;
0540   int index = -1;
0541   for (unsigned int i = 0; i < input.size(); ++i) {
0542     int quality = input[i].getQuality();
0543     if (quality > bestQ) {
0544       bestQ = quality;
0545       index = i;
0546     }
0547   }
0548   if (index != -1)
0549     return input[index];
0550   return CSCCorrelatedLCTDigi();
0551 }
0552
0553 float CSCStubMatcher::zpositionOfLayer(unsigned int detid, int layer) const {
0554   const auto& id = CSCDetId(detid);
0555   const auto& chamber(cscGeometry_->chamber(id));
0556   return fabs(chamber->layer(layer)->centerOfStrip(20).z());
0557 }
0558
0559 int CSCStubMatcher::nChambersWithCLCT(int min_quality) const {
0560   int result = 0;
0561   const auto& chamber_ids = chamberIdsCLCT();
0562   for (const auto& id : chamber_ids) {
0563     int nStubChamber = 0;
0564     const auto& clcts = clctsInChamber(id);
0565     for (const auto& clct : clcts) {
0566       if (!clct.isValid())
0567         continue;
0568       if (clct.getQuality() >= min_quality) {
0569         nStubChamber++;
0570       }
0571     }
0572     if (nStubChamber > 0) {
0573       ++result;
0574     }
0575   }
0576   return result;
0577 }
0578
0579 int CSCStubMatcher::nChambersWithALCT(int min_quality) const {
0580   int result = 0;
0581   const auto& chamber_ids = chamberIdsALCT();
0582   for (const auto& id : chamber_ids) {
0583     int nStubChamber = 0;
0584     const auto& alcts = alctsInChamber(id);
0585     for (const auto& alct : alcts) {
0586       if (!alct.isValid())
0587         continue;
0588       if (alct.getQuality() >= min_quality) {
0589         nStubChamber++;
0590       }
0591     }
0592     if (nStubChamber > 0) {
0593       ++result;
0594     }
0595   }
0596   return result;
0597 }
0598
0599 int CSCStubMatcher::nChambersWithLCT(int min_quality) const {
0600   int result = 0;
0601   const auto& chamber_ids = chamberIdsLCT();
0602   for (const auto& id : chamber_ids) {
0603     int nStubChamber = 0;
0604     const auto& lcts = lctsInChamber(id);
0605     for (const auto& lct : lcts) {
0606       if (!lct.isValid())
0607         continue;
0608       if (lct.getQuality() >= min_quality) {
0609         nStubChamber++;
0610       }
0611     }
0612     if (nStubChamber > 0) {
0613       ++result;
0614     }
0615   }
0616   return result;
0617 }
0618
0619 int CSCStubMatcher::nChambersWithMPLCT(int min_quality) const {
0620   int result = 0;
0621   const auto& chamber_ids = chamberIdsMPLCT();
0622   for (const auto& id : chamber_ids) {
0623     int nStubChamber = 0;
0624     const auto& mplcts = mplctsInChamber(id);
0625     for (const auto& mplct : mplcts) {
0626       if (!mplct.isValid())
0627         continue;
0628       if (mplct.getQuality() >= min_quality) {
0629         nStubChamber++;
0630       }
0631     }
0632     if (nStubChamber > 0) {
0633       ++result;
0634     }
0635   }
0636   return result;
0637 }
0638
0639 bool CSCStubMatcher::lctInChamber(const CSCDetId& id, const CSCCorrelatedLCTDigi& lct) const {
0640   for (const auto& stub : lctsInChamber(id.rawId())) {
0641     if (stub == lct)
0642       return true;
0643   }
0644   return false;
0645 }
0646
0647 GlobalPoint CSCStubMatcher::getGlobalPosition(unsigned int rawId, const CSCCorrelatedLCTDigi& lct) const {
0648   CSCDetId cscId(rawId);
0649   CSCDetId keyId(cscId.endcap(), cscId.station(), cscId.ring(), cscId.chamber(), CSCConstants::KEY_CLCT_LAYER);
0650   float fractional_strip = lct.getFractionalStrip();
0651   // case ME1/1
0652   if (cscId.station() == 1 and (cscId.ring() == 4 || cscId.ring() == 1)) {
0653     int ring = 1;  // Default to ME1/b
0654     if (lct.getStrip() > CSCConstants::MAX_HALF_STRIP_ME1B) {
0655       ring = 4;  // Change to ME1/a if the HalfStrip Number exceeds the range of ME1/b
0656       fractional_strip -= CSCConstants::NUM_STRIPS_ME1B;
0657     }
0658     CSCDetId cscId_(cscId.endcap(), cscId.station(), ring, cscId.chamber(), cscId.layer());
0659     cscId = cscId_;
0660   }
0661   // regular cases
0662   const auto& chamber = cscGeometry_->chamber(cscId);
0663   const auto& layer_geo = chamber->layer(CSCConstants::KEY_CLCT_LAYER)->geometry();
0664   // LCT::getKeyWG() also starts from 0
0665   float wire = layer_geo->middleWireOfGroup(lct.getKeyWG() + 1);
0666   const LocalPoint& csc_intersect = layer_geo->intersectionOfStripAndWire(fractional_strip, wire);
0667   const GlobalPoint& csc_gp = cscGeometry_->idToDet(keyId)->surface().toGlobal(csc_intersect);
0668   return csc_gp;
0669 }
0670
0671 void CSCStubMatcher::clear() {
0672   chamber_to_clcts_all_.clear();
0673   chamber_to_alcts_all_.clear();
0674   chamber_to_lcts_all_.clear();
0675   chamber_to_mplcts_all_.clear();
0676
0677   chamber_to_clcts_.clear();
0678   chamber_to_alcts_.clear();
0679   chamber_to_lcts_.clear();
0680   chamber_to_mplcts_.clear();
0681 }
0682
0683 void CSCStubMatcher::addGhostLCTs(const CSCCorrelatedLCTDigi& lct11,
0684                                   const CSCCorrelatedLCTDigi& lct22,
0685                                   CSCCorrelatedLCTDigiContainer& lcts_tmp) const {
0686   int wg1 = lct11.getKeyWG();
0687   int wg2 = lct22.getKeyWG();
0688   int hs1 = lct11.getStrip();
0689   int hs2 = lct22.getStrip();
0690
0691   if (!(wg1 == wg2 || hs1 == hs2)) {
0692     // flip the ALCTs
0693     CSCCorrelatedLCTDigi lct12 = lct11;
0694     lct12.setWireGroup(wg2);
0695     lct12.setALCT(lct22.getALCT());
0696     lct12.setCLCT(lct11.getCLCT());
0697     lcts_tmp.push_back(lct12);
0698
0699     CSCCorrelatedLCTDigi lct21 = lct22;
0700     lct21.setWireGroup(wg1);
0701     lct21.setALCT(lct11.getALCT());
0702     lct21.setCLCT(lct22.getCLCT());
0703     lcts_tmp.push_back(lct21);
0704   }
0705 }