MuonReco/src/Muon.cc

0001 #include "DataFormats/MuonReco/interface/Muon.h"
0002 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
0003 #include "DataFormats/MuonDetId/interface/RPCDetId.h"
0004 #include "DataFormats/MuonDetId/interface/GEMDetId.h"
0005 #include "DataFormats/MuonDetId/interface/ME0DetId.h"
0006
0007 using namespace reco;
0008
0009 Muon::Muon(Charge q, const LorentzVector& p4, const Point& vtx) : RecoCandidate(q, p4, vtx, -13 * q) {
0010   energyValid_ = false;
0011   matchesValid_ = false;
0012   isolationValid_ = false;
0013   pfIsolationValid_ = false;
0014   qualityValid_ = false;
0015   caloCompatibility_ = -9999.;
0016   type_ = 0;
0017   bestTunePTrackType_ = reco::Muon::None;
0018   bestTrackType_ = reco::Muon::None;
0019   selectors_ = 0;
0020 }
0021
0022 Muon::Muon() {
0023   energyValid_ = false;
0024   matchesValid_ = false;
0025   isolationValid_ = false;
0026   pfIsolationValid_ = false;
0027   qualityValid_ = false;
0028   caloCompatibility_ = -9999.;
0029   type_ = 0;
0030   bestTrackType_ = reco::Muon::None;
0031   bestTunePTrackType_ = reco::Muon::None;
0032   selectors_ = 0;
0033 }
0034
0035 bool Muon::overlap(const Candidate& c) const {
0036   const RecoCandidate* o = dynamic_cast<const RecoCandidate*>(&c);
0037   return (o != nullptr && (checkOverlap(track(), o->track()) || checkOverlap(standAloneMuon(), o->standAloneMuon()) ||
0038                            checkOverlap(combinedMuon(), o->combinedMuon()) ||
0039                            checkOverlap(standAloneMuon(), o->track()) || checkOverlap(combinedMuon(), o->track())));
0040 }
0041
0042 Muon* Muon::clone() const { return new Muon(*this); }
0043
0044 int Muon::numberOfChambersCSCorDT() const {
0045   int total = 0;
0046   int nAll = numberOfChambers();
0047   for (int iC = 0; iC < nAll; ++iC) {
0048     if (matches()[iC].detector() == MuonSubdetId::CSC || matches()[iC].detector() == MuonSubdetId::DT)
0049       total++;
0050   }
0051
0052   return total;
0053 }
0054
0055 int Muon::numberOfMatches(ArbitrationType type) const {
0056   int matches(0);
0057   for (auto& chamberMatch : muMatches_) {
0058     if (type == RPCHitAndTrackArbitration) {
0059       if (chamberMatch.rpcMatches.empty())
0060         continue;
0061       matches += chamberMatch.rpcMatches.size();
0062       continue;
0063     }
0064     if (type == ME0SegmentAndTrackArbitration) {
0065       if (chamberMatch.me0Matches.empty())
0066         continue;
0067       matches += chamberMatch.me0Matches.size();
0068       continue;
0069     }
0070     if (type == GEMSegmentAndTrackArbitration) {
0071       if (chamberMatch.gemMatches.empty())
0072         continue;
0073       matches += chamberMatch.gemMatches.size();
0074       continue;
0075     }
0076
0077     if (chamberMatch.segmentMatches.empty())
0078       continue;
0079     if (type == NoArbitration) {
0080       matches++;
0081       continue;
0082     }
0083
0084     for (auto& segmentMatch : chamberMatch.segmentMatches) {
0085       if (type == SegmentArbitration)
0086         if (segmentMatch.isMask(MuonSegmentMatch::BestInChamberByDR)) {
0087           matches++;
0088           break;
0089         }
0090       if (type == SegmentAndTrackArbitration)
0091         if (segmentMatch.isMask(MuonSegmentMatch::BestInChamberByDR) &&
0092             segmentMatch.isMask(MuonSegmentMatch::BelongsToTrackByDR)) {
0093           matches++;
0094           break;
0095         }
0096       if (type == SegmentAndTrackArbitrationCleaned)
0097         if (segmentMatch.isMask(MuonSegmentMatch::BestInChamberByDR) &&
0098             segmentMatch.isMask(MuonSegmentMatch::BelongsToTrackByDR) &&
0099             segmentMatch.isMask(MuonSegmentMatch::BelongsToTrackByCleaning)) {
0100           matches++;
0101           break;
0102         }
0103     }
0104   }
0105
0106   return matches;
0107 }
0108
0109 int Muon::numberOfMatchedStations(ArbitrationType type) const {
0110   int stations(0);
0111
0112   unsigned int theStationMask = stationMask(type);
0113   // eight stations, eight bits
0114   for (int it = 0; it < 8; ++it)
0115     if (theStationMask & 1 << it)
0116       ++stations;
0117
0118   return stations;
0119 }
0120
0121 unsigned int Muon::expectedNnumberOfMatchedStations(float minDistanceFromEdge) const {
0122   unsigned int stationMask = 0;
0123   minDistanceFromEdge = std::abs(minDistanceFromEdge);
0124   for (auto& chamberMatch : muMatches_) {
0125     if (chamberMatch.detector() != MuonSubdetId::DT && chamberMatch.detector() != MuonSubdetId::CSC)
0126       continue;
0127     float edgeX = chamberMatch.edgeX;
0128     float edgeY = chamberMatch.edgeY;
0129     // check we if the trajectory is well within the acceptance
0130     if (edgeX < 0 && -edgeX > minDistanceFromEdge && edgeY < 0 && -edgeY > minDistanceFromEdge)
0131       stationMask |= 1 << ((chamberMatch.station() - 1) + 4 * (chamberMatch.detector() - 1));
0132   }
0133   unsigned int n = 0;
0134   for (unsigned int i = 0; i < 8; ++i)
0135     if (stationMask & (1 << i))
0136       n++;
0137   return n;
0138 }
0139
0140 unsigned int Muon::stationMask(ArbitrationType type) const {
0141   unsigned int totMask(0);
0142   unsigned int curMask(0);
0143
0144   for (auto& chamberMatch : muMatches_) {
0145     if (type == RPCHitAndTrackArbitration) {
0146       if (chamberMatch.rpcMatches.empty())
0147         continue;
0148       RPCDetId rollId = chamberMatch.id.rawId();
0149       int rpcIndex = rollId.region() == 0 ? 1 : 2;
0150       curMask = 1 << ((chamberMatch.station() - 1) + 4 * (rpcIndex - 1));
0151       // do not double count
0152       if (!(totMask & curMask))
0153         totMask += curMask;
0154       continue;
0155     }
0156
0157     if (chamberMatch.segmentMatches.empty())
0158       continue;
0159     if (type == NoArbitration) {
0160       curMask = 1 << ((chamberMatch.station() - 1) + 4 * (chamberMatch.detector() - 1));
0161       // do not double count
0162       if (!(totMask & curMask))
0163         totMask += curMask;
0164       continue;
0165     }
0166
0167     for (auto& segmentMatch : chamberMatch.segmentMatches) {
0168       if (type == SegmentArbitration)
0169         if (segmentMatch.isMask(MuonSegmentMatch::BestInStationByDR)) {
0170           curMask = 1 << ((chamberMatch.station() - 1) + 4 * (chamberMatch.detector() - 1));
0171           // do not double count
0172           if (!(totMask & curMask))
0173             totMask += curMask;
0174           break;
0175         }
0176       if (type == SegmentAndTrackArbitration)
0177         if (segmentMatch.isMask(MuonSegmentMatch::BestInStationByDR) &&
0178             segmentMatch.isMask(MuonSegmentMatch::BelongsToTrackByDR)) {
0179           curMask = 1 << ((chamberMatch.station() - 1) + 4 * (chamberMatch.detector() - 1));
0180           // do not double count
0181           if (!(totMask & curMask))
0182             totMask += curMask;
0183           break;
0184         }
0185       if (type == SegmentAndTrackArbitrationCleaned)
0186         if (segmentMatch.isMask(MuonSegmentMatch::BestInStationByDR) &&
0187             segmentMatch.isMask(MuonSegmentMatch::BelongsToTrackByDR) &&
0188             segmentMatch.isMask(MuonSegmentMatch::BelongsToTrackByCleaning)) {
0189           curMask = 1 << ((chamberMatch.station() - 1) + 4 * (chamberMatch.detector() - 1));
0190           // do not double count
0191           if (!(totMask & curMask))
0192             totMask += curMask;
0193           break;
0194         }
0195     }
0196   }
0197
0198   return totMask;
0199 }
0200
0201 int Muon::numberOfMatchedRPCLayers(ArbitrationType type) const {
0202   int layers(0);
0203
0204   unsigned int theRPCLayerMask = RPClayerMask(type);
0205   // maximum ten layers because of 6 layers in barrel and 3 (4) layers in each endcap before (after) upscope
0206   for (int it = 0; it < 10; ++it)
0207     if (theRPCLayerMask & 1 << it)
0208       ++layers;
0209
0210   return layers;
0211 }
0212
0213 unsigned int Muon::RPClayerMask(ArbitrationType type) const {
0214   unsigned int totMask(0);
0215   unsigned int curMask(0);
0216   for (auto& chamberMatch : muMatches_) {
0217     if (chamberMatch.rpcMatches.empty())
0218       continue;
0219
0220     RPCDetId rollId = chamberMatch.id.rawId();
0221     const int region = rollId.region();
0222     const int stationIndex = chamberMatch.station();
0223     int rpcLayer = stationIndex;
0224     if (region == 0) {
0225       const int layer = rollId.layer();
0226       rpcLayer = stationIndex - 1 + stationIndex * layer;
0227       if ((stationIndex == 2 && layer == 2) || (stationIndex == 4 && layer == 1))
0228         rpcLayer -= 1;
0229     } else
0230       rpcLayer += 6;
0231
0232     curMask = 1 << (rpcLayer - 1);
0233     // do not double count
0234     if (!(totMask & curMask))
0235       totMask += curMask;
0236   }
0237
0238   return totMask;
0239 }
0240
0241 unsigned int Muon::stationGapMaskDistance(float distanceCut) const {
0242   unsigned int totMask(0);
0243   distanceCut = std::abs(distanceCut);
0244   for (auto& chamberMatch : muMatches_) {
0245     unsigned int curMask(0);
0246     const int detectorIndex = chamberMatch.detector();
0247     if (detectorIndex < 1 || detectorIndex >= 4)
0248       continue;
0249     const int stationIndex = chamberMatch.station();
0250     if (stationIndex < 1 || stationIndex >= 5)
0251       continue;
0252
0253     float edgeX = chamberMatch.edgeX;
0254     float edgeY = chamberMatch.edgeY;
0255     if (edgeX < 0 && -edgeX > distanceCut && edgeY < 0 &&
0256         -edgeY > distanceCut)  // inside the chamber so negates all gaps for this station
0257       continue;
0258
0259     if ((std::abs(edgeX) < distanceCut && edgeY < distanceCut) ||
0260         (std::abs(edgeY) < distanceCut && edgeX < distanceCut))  // inside gap
0261       curMask = 1 << ((stationIndex - 1) + 4 * (detectorIndex - 1));
0262     totMask += curMask;  // add to total mask
0263   }
0264
0265   return totMask;
0266 }
0267
0268 unsigned int Muon::stationGapMaskPull(float sigmaCut) const {
0269   sigmaCut = std::abs(sigmaCut);
0270   unsigned int totMask(0);
0271   for (auto& chamberMatch : muMatches_) {
0272     unsigned int curMask(0);
0273     const int detectorIndex = chamberMatch.detector();
0274     if (detectorIndex < 1 || detectorIndex >= 4)
0275       continue;
0276     const int stationIndex = chamberMatch.station();
0277     if (stationIndex < 1 || stationIndex >= 5)
0278       continue;
0279
0280     float edgeX = chamberMatch.edgeX;
0281     float edgeY = chamberMatch.edgeY;
0282     float xErr = chamberMatch.xErr + 0.000001;  // protect against division by zero later
0283     float yErr = chamberMatch.yErr + 0.000001;  // protect against division by zero later
0284     if (edgeX < 0 && std::abs(edgeX / xErr) > sigmaCut && edgeY < 0 &&
0285         std::abs(edgeY / yErr) > sigmaCut)  // inside the chamber so negates all gaps for this station
0286       continue;
0287
0288     if ((std::abs(edgeX / xErr) < sigmaCut && edgeY < sigmaCut * yErr) ||
0289         (std::abs(edgeY / yErr) < sigmaCut && edgeX < sigmaCut * xErr))  // inside gap
0290       curMask = 1 << ((stationIndex - 1) + 4 * (detectorIndex - 1));
0291     totMask += curMask;  // add to total mask
0292   }
0293
0294   return totMask;
0295 }
0296
0297 int Muon::nDigisInStation(int station, int muonSubdetId) const {
0298   int nDigis(0);
0299   std::map<int, int> me11DigisPerCh;
0300
0301   if (muonSubdetId != MuonSubdetId::CSC && muonSubdetId != MuonSubdetId::DT)
0302     return 0;
0303
0304   for (auto& match : muMatches_) {
0305     if (match.detector() != muonSubdetId || match.station() != station)
0306       continue;
0307
0308     int nDigisInCh = match.nDigisInRange;
0309
0310     if (muonSubdetId == MuonSubdetId::CSC && station == 1) {
0311       CSCDetId id(match.id.rawId());
0312
0313       int chamber = id.chamber();
0314       int ring = id.ring();
0315
0316       if (ring == 1 || ring == 4)  // merge ME1/1a and ME1/1b digis
0317       {
0318         if (me11DigisPerCh.find(chamber) == me11DigisPerCh.end())
0319           me11DigisPerCh[chamber] = 0;
0320
0321         me11DigisPerCh[chamber] += nDigisInCh;
0322
0323         continue;
0324       }
0325     }
0326
0327     if (nDigisInCh > nDigis)
0328       nDigis = nDigisInCh;
0329   }
0330
0331   for (const auto& me11DigisInCh : me11DigisPerCh) {
0332     int nMe11DigisInCh = me11DigisInCh.second;
0333     if (nMe11DigisInCh > nDigis)
0334       nDigis = nMe11DigisInCh;
0335   }
0336
0337   return nDigis;
0338 }
0339
0340 bool Muon::hasShowerInStation(int station, int muonSubdetId, int nDtDigisCut, int nCscDigisCut) const {
0341   if (muonSubdetId != MuonSubdetId::DT && muonSubdetId != MuonSubdetId::CSC)
0342     return false;
0343   auto nDigisCut = muonSubdetId == MuonSubdetId::DT ? nDtDigisCut : nCscDigisCut;
0344
0345   return nDigisInStation(station, muonSubdetId) >= nDigisCut;
0346 }
0347
0348 int Muon::numberOfShowers(int nDtDigisCut, int nCscDigisCut) const {
0349   int nShowers = 0;
0350   for (int station = 1; station < 5; ++station) {
0351     if (hasShowerInStation(station, MuonSubdetId::DT, nDtDigisCut, nCscDigisCut))
0352       nShowers++;
0353     if (hasShowerInStation(station, MuonSubdetId::CSC, nDtDigisCut, nCscDigisCut))
0354       nShowers++;
0355   }
0356
0357   return nShowers;
0358 }
0359
0360 int Muon::numberOfSegments(int station, int muonSubdetId, ArbitrationType type) const {
0361   int segments(0);
0362   for (auto& chamberMatch : muMatches_) {
0363     if (chamberMatch.segmentMatches.empty())
0364       continue;
0365     if (chamberMatch.detector() != muonSubdetId || chamberMatch.station() != station)
0366       continue;
0367
0368     if (type == NoArbitration) {
0369       segments += chamberMatch.segmentMatches.size();
0370       continue;
0371     }
0372
0373     for (auto& segmentMatch : chamberMatch.segmentMatches) {
0374       if (type == SegmentArbitration)
0375         if (segmentMatch.isMask(MuonSegmentMatch::BestInStationByDR)) {
0376           segments++;
0377           break;
0378         }
0379       if (type == SegmentAndTrackArbitration)
0380         if (segmentMatch.isMask(MuonSegmentMatch::BestInStationByDR) &&
0381             segmentMatch.isMask(MuonSegmentMatch::BelongsToTrackByDR)) {
0382           segments++;
0383           break;
0384         }
0385       if (type == SegmentAndTrackArbitrationCleaned)
0386         if (segmentMatch.isMask(MuonSegmentMatch::BestInStationByDR) &&
0387             segmentMatch.isMask(MuonSegmentMatch::BelongsToTrackByDR) &&
0388             segmentMatch.isMask(MuonSegmentMatch::BelongsToTrackByCleaning)) {
0389           segments++;
0390           break;
0391         }
0392     }
0393   }
0394
0395   return segments;
0396 }
0397
0398 const std::vector<const MuonChamberMatch*> Muon::chambers(int station, int muonSubdetId) const {
0399   std::vector<const MuonChamberMatch*> chambers;
0400   for (std::vector<MuonChamberMatch>::const_iterator chamberMatch = muMatches_.begin();
0401        chamberMatch != muMatches_.end();
0402        chamberMatch++)
0403     if (chamberMatch->detector() == muonSubdetId && chamberMatch->station() == station)
0404       chambers.push_back(&(*chamberMatch));
0405   return chambers;
0406 }
0407
0408 std::pair<const MuonChamberMatch*, const MuonSegmentMatch*> Muon::pair(
0409     const std::vector<const MuonChamberMatch*>& chambers, ArbitrationType type) const {
0410   MuonChamberMatch* m = nullptr;
0411   MuonSegmentMatch* s = nullptr;
0412   std::pair<const MuonChamberMatch*, const MuonSegmentMatch*> chamberSegmentPair(m, s);
0413
0414   if (chambers.empty())
0415     return chamberSegmentPair;
0416   for (std::vector<const MuonChamberMatch*>::const_iterator chamberMatch = chambers.begin();
0417        chamberMatch != chambers.end();
0418        chamberMatch++) {
0419     if ((*chamberMatch)->segmentMatches.empty())
0420       continue;
0421     if (type == NoArbitration)
0422       return std::make_pair(*chamberMatch, &((*chamberMatch)->segmentMatches.front()));
0423
0424     for (std::vector<MuonSegmentMatch>::const_iterator segmentMatch = (*chamberMatch)->segmentMatches.begin();
0425          segmentMatch != (*chamberMatch)->segmentMatches.end();
0426          segmentMatch++) {
0427       if (type == SegmentArbitration)
0428         if (segmentMatch->isMask(MuonSegmentMatch::BestInStationByDR))
0429           return std::make_pair(*chamberMatch, &(*segmentMatch));
0430       if (type == SegmentAndTrackArbitration)
0431         if (segmentMatch->isMask(MuonSegmentMatch::BestInStationByDR) &&
0432             segmentMatch->isMask(MuonSegmentMatch::BelongsToTrackByDR))
0433           return std::make_pair(*chamberMatch, &(*segmentMatch));
0434       if (type == SegmentAndTrackArbitrationCleaned)
0435         if (segmentMatch->isMask(MuonSegmentMatch::BestInStationByDR) &&
0436             segmentMatch->isMask(MuonSegmentMatch::BelongsToTrackByDR) &&
0437             segmentMatch->isMask(MuonSegmentMatch::BelongsToTrackByCleaning))
0438           return std::make_pair(*chamberMatch, &(*segmentMatch));
0439     }
0440   }
0441
0442   return chamberSegmentPair;
0443 }
0444
0445 float Muon::dX(int station, int muonSubdetId, ArbitrationType type) const {
0446   std::pair<const MuonChamberMatch*, const MuonSegmentMatch*> chamberSegmentPair =
0447       pair(chambers(station, muonSubdetId), type);
0448   if (chamberSegmentPair.first == nullptr || chamberSegmentPair.second == nullptr)
0449     return 999999;
0450   if (!chamberSegmentPair.second->hasPhi())
0451     return 999999;
0452   return chamberSegmentPair.first->x - chamberSegmentPair.second->x;
0453 }
0454
0455 float Muon::dY(int station, int muonSubdetId, ArbitrationType type) const {
0456   if (station == 4 && muonSubdetId == MuonSubdetId::DT)
0457     return 999999;  // no y information
0458   std::pair<const MuonChamberMatch*, const MuonSegmentMatch*> chamberSegmentPair =
0459       pair(chambers(station, muonSubdetId), type);
0460   if (chamberSegmentPair.first == nullptr || chamberSegmentPair.second == nullptr)
0461     return 999999;
0462   if (!chamberSegmentPair.second->hasZed())
0463     return 999999;
0464   return chamberSegmentPair.first->y - chamberSegmentPair.second->y;
0465 }
0466
0467 float Muon::dDxDz(int station, int muonSubdetId, ArbitrationType type) const {
0468   std::pair<const MuonChamberMatch*, const MuonSegmentMatch*> chamberSegmentPair =
0469       pair(chambers(station, muonSubdetId), type);
0470   if (chamberSegmentPair.first == nullptr || chamberSegmentPair.second == nullptr)
0471     return 999999;
0472   if (!chamberSegmentPair.second->hasPhi())
0473     return 999999;
0474   return chamberSegmentPair.first->dXdZ - chamberSegmentPair.second->dXdZ;
0475 }
0476
0477 float Muon::dDyDz(int station, int muonSubdetId, ArbitrationType type) const {
0478   if (station == 4 && muonSubdetId == MuonSubdetId::DT)
0479     return 999999;  // no y information
0480   std::pair<const MuonChamberMatch*, const MuonSegmentMatch*> chamberSegmentPair =
0481       pair(chambers(station, muonSubdetId), type);
0482   if (chamberSegmentPair.first == nullptr || chamberSegmentPair.second == nullptr)
0483     return 999999;
0484   if (!chamberSegmentPair.second->hasZed())
0485     return 999999;
0486   return chamberSegmentPair.first->dYdZ - chamberSegmentPair.second->dYdZ;
0487 }
0488
0489 float Muon::pullX(int station, int muonSubdetId, ArbitrationType type, bool includeSegmentError) const {
0490   std::pair<const MuonChamberMatch*, const MuonSegmentMatch*> chamberSegmentPair =
0491       pair(chambers(station, muonSubdetId), type);
0492   if (chamberSegmentPair.first == nullptr || chamberSegmentPair.second == nullptr)
0493     return 999999;
0494   if (!chamberSegmentPair.second->hasPhi())
0495     return 999999;
0496   if (includeSegmentError)
0497     return (chamberSegmentPair.first->x - chamberSegmentPair.second->x) /
0498            sqrt(std::pow(chamberSegmentPair.first->xErr, 2) + std::pow(chamberSegmentPair.second->xErr, 2));
0499   return (chamberSegmentPair.first->x - chamberSegmentPair.second->x) / chamberSegmentPair.first->xErr;
0500 }
0501
0502 float Muon::pullY(int station, int muonSubdetId, ArbitrationType type, bool includeSegmentError) const {
0503   if (station == 4 && muonSubdetId == MuonSubdetId::DT)
0504     return 999999;  // no y information
0505   std::pair<const MuonChamberMatch*, const MuonSegmentMatch*> chamberSegmentPair =
0506       pair(chambers(station, muonSubdetId), type);
0507   if (chamberSegmentPair.first == nullptr || chamberSegmentPair.second == nullptr)
0508     return 999999;
0509   if (!chamberSegmentPair.second->hasZed())
0510     return 999999;
0511   if (includeSegmentError)
0512     return (chamberSegmentPair.first->y - chamberSegmentPair.second->y) /
0513            sqrt(std::pow(chamberSegmentPair.first->yErr, 2) + std::pow(chamberSegmentPair.second->yErr, 2));
0514   return (chamberSegmentPair.first->y - chamberSegmentPair.second->y) / chamberSegmentPair.first->yErr;
0515 }
0516
0517 float Muon::pullDxDz(int station, int muonSubdetId, ArbitrationType type, bool includeSegmentError) const {
0518   std::pair<const MuonChamberMatch*, const MuonSegmentMatch*> chamberSegmentPair =
0519       pair(chambers(station, muonSubdetId), type);
0520   if (chamberSegmentPair.first == nullptr || chamberSegmentPair.second == nullptr)
0521     return 999999;
0522   if (!chamberSegmentPair.second->hasPhi())
0523     return 999999;
0524   if (includeSegmentError)
0525     return (chamberSegmentPair.first->dXdZ - chamberSegmentPair.second->dXdZ) /
0526            sqrt(std::pow(chamberSegmentPair.first->dXdZErr, 2) + std::pow(chamberSegmentPair.second->dXdZErr, 2));
0527   return (chamberSegmentPair.first->dXdZ - chamberSegmentPair.second->dXdZ) / chamberSegmentPair.first->dXdZErr;
0528 }
0529
0530 float Muon::pullDyDz(int station, int muonSubdetId, ArbitrationType type, bool includeSegmentError) const {
0531   if (station == 4 && muonSubdetId == MuonSubdetId::DT)
0532     return 999999;  // no y information
0533   std::pair<const MuonChamberMatch*, const MuonSegmentMatch*> chamberSegmentPair =
0534       pair(chambers(station, muonSubdetId), type);
0535   if (chamberSegmentPair.first == nullptr || chamberSegmentPair.second == nullptr)
0536     return 999999;
0537   if (!chamberSegmentPair.second->hasZed())
0538     return 999999;
0539   if (includeSegmentError)
0540     return (chamberSegmentPair.first->dYdZ - chamberSegmentPair.second->dYdZ) /
0541            sqrt(std::pow(chamberSegmentPair.first->dYdZErr, 2) + std::pow(chamberSegmentPair.second->dYdZErr, 2));
0542   return (chamberSegmentPair.first->dYdZ - chamberSegmentPair.second->dYdZ) / chamberSegmentPair.first->dYdZErr;
0543 }
0544
0545 float Muon::segmentX(int station, int muonSubdetId, ArbitrationType type) const {
0546   std::pair<const MuonChamberMatch*, const MuonSegmentMatch*> chamberSegmentPair =
0547       pair(chambers(station, muonSubdetId), type);
0548   if (chamberSegmentPair.first == nullptr || chamberSegmentPair.second == nullptr)
0549     return 999999;
0550   if (!chamberSegmentPair.second->hasPhi())
0551     return 999999;
0552   return chamberSegmentPair.second->x;
0553 }
0554
0555 float Muon::segmentY(int station, int muonSubdetId, ArbitrationType type) const {
0556   if (station == 4 && muonSubdetId == MuonSubdetId::DT)
0557     return 999999;  // no y information
0558   std::pair<const MuonChamberMatch*, const MuonSegmentMatch*> chamberSegmentPair =
0559       pair(chambers(station, muonSubdetId), type);
0560   if (chamberSegmentPair.first == nullptr || chamberSegmentPair.second == nullptr)
0561     return 999999;
0562   if (!chamberSegmentPair.second->hasZed())
0563     return 999999;
0564   return chamberSegmentPair.second->y;
0565 }
0566
0567 float Muon::segmentDxDz(int station, int muonSubdetId, ArbitrationType type) const {
0568   std::pair<const MuonChamberMatch*, const MuonSegmentMatch*> chamberSegmentPair =
0569       pair(chambers(station, muonSubdetId), type);
0570   if (chamberSegmentPair.first == nullptr || chamberSegmentPair.second == nullptr)
0571     return 999999;
0572   if (!chamberSegmentPair.second->hasPhi())
0573     return 999999;
0574   return chamberSegmentPair.second->dXdZ;
0575 }
0576
0577 float Muon::segmentDyDz(int station, int muonSubdetId, ArbitrationType type) const {
0578   if (station == 4 && muonSubdetId == MuonSubdetId::DT)
0579     return 999999;  // no y information
0580   std::pair<const MuonChamberMatch*, const MuonSegmentMatch*> chamberSegmentPair =
0581       pair(chambers(station, muonSubdetId), type);
0582   if (chamberSegmentPair.first == nullptr || chamberSegmentPair.second == nullptr)
0583     return 999999;
0584   if (!chamberSegmentPair.second->hasZed())
0585     return 999999;
0586   return chamberSegmentPair.second->dYdZ;
0587 }
0588
0589 float Muon::segmentXErr(int station, int muonSubdetId, ArbitrationType type) const {
0590   std::pair<const MuonChamberMatch*, const MuonSegmentMatch*> chamberSegmentPair =
0591       pair(chambers(station, muonSubdetId), type);
0592   if (chamberSegmentPair.first == nullptr || chamberSegmentPair.second == nullptr)
0593     return 999999;
0594   if (!chamberSegmentPair.second->hasPhi())
0595     return 999999;
0596   return chamberSegmentPair.second->xErr;
0597 }
0598
0599 float Muon::segmentYErr(int station, int muonSubdetId, ArbitrationType type) const {
0600   if (station == 4 && muonSubdetId == MuonSubdetId::DT)
0601     return 999999;  // no y information
0602   std::pair<const MuonChamberMatch*, const MuonSegmentMatch*> chamberSegmentPair =
0603       pair(chambers(station, muonSubdetId), type);
0604   if (chamberSegmentPair.first == nullptr || chamberSegmentPair.second == nullptr)
0605     return 999999;
0606   if (!chamberSegmentPair.second->hasZed())
0607     return 999999;
0608   return chamberSegmentPair.second->yErr;
0609 }
0610
0611 float Muon::segmentDxDzErr(int station, int muonSubdetId, ArbitrationType type) const {
0612   std::pair<const MuonChamberMatch*, const MuonSegmentMatch*> chamberSegmentPair =
0613       pair(chambers(station, muonSubdetId), type);
0614   if (chamberSegmentPair.first == nullptr || chamberSegmentPair.second == nullptr)
0615     return 999999;
0616   if (!chamberSegmentPair.second->hasPhi())
0617     return 999999;
0618   return chamberSegmentPair.second->dXdZErr;
0619 }
0620
0621 float Muon::segmentDyDzErr(int station, int muonSubdetId, ArbitrationType type) const {
0622   if (station == 4 && muonSubdetId == MuonSubdetId::DT)
0623     return 999999;  // no y information
0624   std::pair<const MuonChamberMatch*, const MuonSegmentMatch*> chamberSegmentPair =
0625       pair(chambers(station, muonSubdetId), type);
0626   if (chamberSegmentPair.first == nullptr || chamberSegmentPair.second == nullptr)
0627     return 999999;
0628   if (!chamberSegmentPair.second->hasZed())
0629     return 999999;
0630   return chamberSegmentPair.second->dYdZErr;
0631 }
0632
0633 float Muon::trackEdgeX(int station, int muonSubdetId, ArbitrationType type) const {
0634   const std::vector<const MuonChamberMatch*> muonChambers = chambers(station, muonSubdetId);
0635   if (muonChambers.empty())
0636     return 999999;
0637
0638   std::pair<const MuonChamberMatch*, const MuonSegmentMatch*> chamberSegmentPair = pair(muonChambers, type);
0639   if (chamberSegmentPair.first == nullptr || chamberSegmentPair.second == nullptr) {
0640     float dist = 999999;
0641     float supVar = 999999;
0642     for (const MuonChamberMatch* muonChamber : muonChambers) {
0643       float currDist = muonChamber->dist();
0644       if (currDist < dist) {
0645         dist = currDist;
0646         supVar = muonChamber->edgeX;
0647       }
0648     }
0649     return supVar;
0650   } else
0651     return chamberSegmentPair.first->edgeX;
0652 }
0653
0654 float Muon::trackEdgeY(int station, int muonSubdetId, ArbitrationType type) const {
0655   const std::vector<const MuonChamberMatch*> muonChambers = chambers(station, muonSubdetId);
0656   if (muonChambers.empty())
0657     return 999999;
0658
0659   std::pair<const MuonChamberMatch*, const MuonSegmentMatch*> chamberSegmentPair = pair(muonChambers, type);
0660   if (chamberSegmentPair.first == nullptr || chamberSegmentPair.second == nullptr) {
0661     float dist = 999999;
0662     float supVar = 999999;
0663     for (const MuonChamberMatch* muonChamber : muonChambers) {
0664       float currDist = muonChamber->dist();
0665       if (currDist < dist) {
0666         dist = currDist;
0667         supVar = muonChamber->edgeY;
0668       }
0669     }
0670     return supVar;
0671   } else
0672     return chamberSegmentPair.first->edgeY;
0673 }
0674
0675 float Muon::trackX(int station, int muonSubdetId, ArbitrationType type) const {
0676   const std::vector<const MuonChamberMatch*> muonChambers = chambers(station, muonSubdetId);
0677   if (muonChambers.empty())
0678     return 999999;
0679
0680   std::pair<const MuonChamberMatch*, const MuonSegmentMatch*> chamberSegmentPair = pair(muonChambers, type);
0681   if (chamberSegmentPair.first == nullptr || chamberSegmentPair.second == nullptr) {
0682     float dist = 999999;
0683     float supVar = 999999;
0684     for (const MuonChamberMatch* muonChamber : muonChambers) {
0685       float currDist = muonChamber->dist();
0686       if (currDist < dist) {
0687         dist = currDist;
0688         supVar = muonChamber->x;
0689       }
0690     }
0691     return supVar;
0692   } else
0693     return chamberSegmentPair.first->x;
0694 }
0695
0696 float Muon::trackY(int station, int muonSubdetId, ArbitrationType type) const {
0697   const std::vector<const MuonChamberMatch*> muonChambers = chambers(station, muonSubdetId);
0698   if (muonChambers.empty())
0699     return 999999;
0700
0701   std::pair<const MuonChamberMatch*, const MuonSegmentMatch*> chamberSegmentPair = pair(muonChambers, type);
0702   if (chamberSegmentPair.first == nullptr || chamberSegmentPair.second == nullptr) {
0703     float dist = 999999;
0704     float supVar = 999999;
0705     for (const MuonChamberMatch* muonChamber : muonChambers) {
0706       float currDist = muonChamber->dist();
0707       if (currDist < dist) {
0708         dist = currDist;
0709         supVar = muonChamber->y;
0710       }
0711     }
0712     return supVar;
0713   } else
0714     return chamberSegmentPair.first->y;
0715 }
0716
0717 float Muon::trackDxDz(int station, int muonSubdetId, ArbitrationType type) const {
0718   const std::vector<const MuonChamberMatch*> muonChambers = chambers(station, muonSubdetId);
0719   if (muonChambers.empty())
0720     return 999999;
0721
0722   std::pair<const MuonChamberMatch*, const MuonSegmentMatch*> chamberSegmentPair = pair(muonChambers, type);
0723   if (chamberSegmentPair.first == nullptr || chamberSegmentPair.second == nullptr) {
0724     float dist = 999999;
0725     float supVar = 999999;
0726     for (const MuonChamberMatch* muonChamber : muonChambers) {
0727       float currDist = muonChamber->dist();
0728       if (currDist < dist) {
0729         dist = currDist;
0730         supVar = muonChamber->dXdZ;
0731       }
0732     }
0733     return supVar;
0734   } else
0735     return chamberSegmentPair.first->dXdZ;
0736 }
0737
0738 float Muon::trackDyDz(int station, int muonSubdetId, ArbitrationType type) const {
0739   const std::vector<const MuonChamberMatch*> muonChambers = chambers(station, muonSubdetId);
0740   if (muonChambers.empty())
0741     return 999999;
0742
0743   std::pair<const MuonChamberMatch*, const MuonSegmentMatch*> chamberSegmentPair = pair(muonChambers, type);
0744   if (chamberSegmentPair.first == nullptr || chamberSegmentPair.second == nullptr) {
0745     float dist = 999999;
0746     float supVar = 999999;
0747     for (const MuonChamberMatch* muonChamber : muonChambers) {
0748       float currDist = muonChamber->dist();
0749       if (currDist < dist) {
0750         dist = currDist;
0751         supVar = muonChamber->dYdZ;
0752       }
0753     }
0754     return supVar;
0755   } else
0756     return chamberSegmentPair.first->dYdZ;
0757 }
0758
0759 float Muon::trackXErr(int station, int muonSubdetId, ArbitrationType type) const {
0760   const std::vector<const MuonChamberMatch*> muonChambers = chambers(station, muonSubdetId);
0761   if (muonChambers.empty())
0762     return 999999;
0763
0764   std::pair<const MuonChamberMatch*, const MuonSegmentMatch*> chamberSegmentPair = pair(muonChambers, type);
0765   if (chamberSegmentPair.first == nullptr || chamberSegmentPair.second == nullptr) {
0766     float dist = 999999;
0767     float supVar = 999999;
0768     for (const MuonChamberMatch* muonChamber : muonChambers) {
0769       float currDist = muonChamber->dist();
0770       if (currDist < dist) {
0771         dist = currDist;
0772         supVar = muonChamber->xErr;
0773       }
0774     }
0775     return supVar;
0776   } else
0777     return chamberSegmentPair.first->xErr;
0778 }
0779
0780 float Muon::trackYErr(int station, int muonSubdetId, ArbitrationType type) const {
0781   const std::vector<const MuonChamberMatch*> muonChambers = chambers(station, muonSubdetId);
0782   if (muonChambers.empty())
0783     return 999999;
0784
0785   std::pair<const MuonChamberMatch*, const MuonSegmentMatch*> chamberSegmentPair = pair(muonChambers, type);
0786   if (chamberSegmentPair.first == nullptr || chamberSegmentPair.second == nullptr) {
0787     float dist = 999999;
0788     float supVar = 999999;
0789     for (const MuonChamberMatch* muonChamber : muonChambers) {
0790       float currDist = muonChamber->dist();
0791       if (currDist < dist) {
0792         dist = currDist;
0793         supVar = muonChamber->yErr;
0794       }
0795     }
0796     return supVar;
0797   } else
0798     return chamberSegmentPair.first->yErr;
0799 }
0800
0801 float Muon::trackDxDzErr(int station, int muonSubdetId, ArbitrationType type) const {
0802   const std::vector<const MuonChamberMatch*> muonChambers = chambers(station, muonSubdetId);
0803   if (muonChambers.empty())
0804     return 999999;
0805
0806   std::pair<const MuonChamberMatch*, const MuonSegmentMatch*> chamberSegmentPair = pair(muonChambers, type);
0807   if (chamberSegmentPair.first == nullptr || chamberSegmentPair.second == nullptr) {
0808     float dist = 999999;
0809     float supVar = 999999;
0810     for (const MuonChamberMatch* muonChamber : muonChambers) {
0811       float currDist = muonChamber->dist();
0812       if (currDist < dist) {
0813         dist = currDist;
0814         supVar = muonChamber->dXdZErr;
0815       }
0816     }
0817     return supVar;
0818   } else
0819     return chamberSegmentPair.first->dXdZErr;
0820 }
0821
0822 float Muon::trackDyDzErr(int station, int muonSubdetId, ArbitrationType type) const {
0823   const std::vector<const MuonChamberMatch*> muonChambers = chambers(station, muonSubdetId);
0824   if (muonChambers.empty())
0825     return 999999;
0826
0827   std::pair<const MuonChamberMatch*, const MuonSegmentMatch*> chamberSegmentPair = pair(muonChambers, type);
0828   if (chamberSegmentPair.first == nullptr || chamberSegmentPair.second == nullptr) {
0829     float dist = 999999;
0830     float supVar = 999999;
0831     for (const MuonChamberMatch* muonChamber : muonChambers) {
0832       float currDist = muonChamber->dist();
0833       if (currDist < dist) {
0834         dist = currDist;
0835         supVar = muonChamber->dYdZErr;
0836       }
0837     }
0838     return supVar;
0839   } else
0840     return chamberSegmentPair.first->dYdZErr;
0841 }
0842
0843 float Muon::trackDist(int station, int muonSubdetId, ArbitrationType type) const {
0844   const std::vector<const MuonChamberMatch*> muonChambers = chambers(station, muonSubdetId);
0845   if (muonChambers.empty())
0846     return 999999;
0847
0848   std::pair<const MuonChamberMatch*, const MuonSegmentMatch*> chamberSegmentPair = pair(muonChambers, type);
0849   if (chamberSegmentPair.first == nullptr || chamberSegmentPair.second == nullptr) {
0850     float dist = 999999;
0851     for (const MuonChamberMatch* muonChamber : muonChambers) {
0852       float currDist = muonChamber->dist();
0853       if (currDist < dist)
0854         dist = currDist;
0855     }
0856     return dist;
0857   } else
0858     return chamberSegmentPair.first->dist();
0859 }
0860
0861 float Muon::trackDistErr(int station, int muonSubdetId, ArbitrationType type) const {
0862   const std::vector<const MuonChamberMatch*> muonChambers = chambers(station, muonSubdetId);
0863   if (muonChambers.empty())
0864     return 999999;
0865
0866   std::pair<const MuonChamberMatch*, const MuonSegmentMatch*> chamberSegmentPair = pair(muonChambers, type);
0867   if (chamberSegmentPair.first == nullptr || chamberSegmentPair.second == nullptr) {
0868     float dist = 999999;
0869     float supVar = 999999;
0870     for (const MuonChamberMatch* muonChamber : muonChambers) {
0871       float currDist = muonChamber->dist();
0872       if (currDist < dist) {
0873         dist = currDist;
0874         supVar = muonChamber->distErr();
0875       }
0876     }
0877     return supVar;
0878   } else
0879     return chamberSegmentPair.first->distErr();
0880 }
0881
0882 void Muon::setIsolation(const MuonIsolation& isoR03, const MuonIsolation& isoR05) {
0883   isolationR03_ = isoR03;
0884   isolationR05_ = isoR05;
0885   isolationValid_ = true;
0886 }
0887
0888 void Muon::setPFIsolation(const std::string& label, const MuonPFIsolation& deposit) {
0889   if (label == "pfIsolationR03")
0890     pfIsolationR03_ = deposit;
0891
0892   if (label == "pfIsolationR04")
0893     pfIsolationR04_ = deposit;
0894
0895   if (label == "pfIsoMeanDRProfileR03")
0896     pfIsoMeanDRR03_ = deposit;
0897
0898   if (label == "pfIsoMeanDRProfileR04")
0899     pfIsoMeanDRR04_ = deposit;
0900
0901   if (label == "pfIsoSumDRProfileR03")
0902     pfIsoSumDRR03_ = deposit;
0903
0904   if (label == "pfIsoSumDRProfileR04")
0905     pfIsoSumDRR04_ = deposit;
0906
0907   pfIsolationValid_ = true;
0908 }
0909
0910 void Muon::setPFP4(const reco::Candidate::LorentzVector& p4) {
0911   pfP4_ = p4;
0912   type_ = type_ | PFMuon;
0913 }
0914
0915 void Muon::setOuterTrack(const TrackRef& t) { outerTrack_ = t; }
0916 void Muon::setInnerTrack(const TrackRef& t) { innerTrack_ = t; }
0917 void Muon::setTrack(const TrackRef& t) { setInnerTrack(t); }
0918 void Muon::setStandAlone(const TrackRef& t) { setOuterTrack(t); }
0919 void Muon::setGlobalTrack(const TrackRef& t) { globalTrack_ = t; }
0920 void Muon::setCombined(const TrackRef& t) { setGlobalTrack(t); }
0921
0922 bool Muon::isAValidMuonTrack(const MuonTrackType& type) const { return muonTrack(type).isNonnull(); }
0923
0924 TrackRef Muon::muonTrack(const MuonTrackType& type) const {
0925   switch (type) {
0926     case InnerTrack:
0927       return innerTrack();
0928     case OuterTrack:
0929       return standAloneMuon();
0930     case CombinedTrack:
0931       return globalTrack();
0932     case TPFMS:
0933       return tpfmsTrack();
0934     case Picky:
0935       return pickyTrack();
0936     case DYT:
0937       return dytTrack();
0938     default:
0939       return muonTrackFromMap(type);
0940   }
0941 }
0942
0943 void Muon::setMuonTrack(const MuonTrackType& type, const TrackRef& t) {
0944   switch (type) {
0945     case InnerTrack:
0946       setInnerTrack(t);
0947       break;
0948     case OuterTrack:
0949       setStandAlone(t);
0950       break;
0951     case CombinedTrack:
0952       setGlobalTrack(t);
0953       break;
0954     default:
0955       refittedTrackMap_[type] = t;
0956       break;
0957   }
0958 }