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