File indexing completed on 2024-09-08 23:52:00
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012 #include "RecoMuon/MuonIdentification/plugins/MuonIdProducer.h"
0013
0014 #include "DataFormats/Math/interface/deltaPhi.h"
0015 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
0016 #include "DataFormats/MuonDetId/interface/DTChamberId.h"
0017 #include "DataFormats/MuonDetId/interface/CSCDetId.h"
0018 #include "DataFormats/MuonDetId/interface/RPCDetId.h"
0019 #include "DataFormats/MuonDetId/interface/GEMDetId.h"
0020 #include "DataFormats/MuonDetId/interface/ME0DetId.h"
0021 #include "DataFormats/MuonReco/interface/MuonCocktails.h"
0022 #include "DataFormats/MuonReco/interface/MuonTime.h"
0023 #include "DataFormats/MuonReco/interface/MuonTimeExtra.h"
0024 #include "DataFormats/MuonReco/interface/MuonTimeExtraMap.h"
0025 #include "DataFormats/RecoCandidate/interface/IsoDeposit.h"
0026 #include "DataFormats/RecoCandidate/interface/IsoDepositFwd.h"
0027
0028 #include "PhysicsTools/IsolationAlgos/interface/IsoDepositExtractorFactory.h"
0029 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0030
0031 #include "RecoMuon/MuonIdentification/interface/MuonMesh.h"
0032 #include "RecoMuon/MuonIdentification/interface/MuonKinkFinder.h"
0033
0034 MuonIdProducer::MuonIdProducer(const edm::ParameterSet& iConfig)
0035 : geomTokenRun_(esConsumes<edm::Transition::BeginRun>()),
0036 propagatorToken_(esConsumes(edm::ESInputTag("", "SteppingHelixPropagatorAny"))) {
0037 LogTrace("MuonIdentification") << "RecoMuon/MuonIdProducer :: Constructor called";
0038
0039 produces<reco::MuonCollection>();
0040 produces<reco::CaloMuonCollection>();
0041 produces<reco::MuonTimeExtraMap>("combined");
0042 produces<reco::MuonTimeExtraMap>("dt");
0043 produces<reco::MuonTimeExtraMap>("csc");
0044
0045 minPt_ = iConfig.getParameter<double>("minPt");
0046 minP_ = iConfig.getParameter<double>("minP");
0047 minPCaloMuon_ = iConfig.getParameter<double>("minPCaloMuon");
0048 minNumberOfMatches_ = iConfig.getParameter<int>("minNumberOfMatches");
0049 addExtraSoftMuons_ = iConfig.getParameter<bool>("addExtraSoftMuons");
0050 maxAbsEta_ = iConfig.getParameter<double>("maxAbsEta");
0051 maxAbsDx_ = iConfig.getParameter<double>("maxAbsDx");
0052 maxAbsPullX2_ = iConfig.getParameter<double>("maxAbsPullX");
0053 maxAbsPullX2_ *= maxAbsPullX2_;
0054 maxAbsDy_ = iConfig.getParameter<double>("maxAbsDy");
0055 maxAbsPullY2_ = iConfig.getParameter<double>("maxAbsPullY");
0056 maxAbsPullY2_ *= maxAbsPullY2_;
0057 fillCaloCompatibility_ = iConfig.getParameter<bool>("fillCaloCompatibility");
0058 fillEnergy_ = iConfig.getParameter<bool>("fillEnergy");
0059 storeCrossedHcalRecHits_ = iConfig.getParameter<bool>("storeCrossedHcalRecHits");
0060 fillMatching_ = iConfig.getParameter<bool>("fillMatching");
0061 fillIsolation_ = iConfig.getParameter<bool>("fillIsolation");
0062 fillShowerDigis_ = iConfig.getParameter<bool>("fillShowerDigis");
0063 writeIsoDeposits_ = iConfig.getParameter<bool>("writeIsoDeposits");
0064 fillGlobalTrackQuality_ = iConfig.getParameter<bool>("fillGlobalTrackQuality");
0065 fillGlobalTrackRefits_ = iConfig.getParameter<bool>("fillGlobalTrackRefits");
0066 arbitrateTrackerMuons_ = iConfig.getParameter<bool>("arbitrateTrackerMuons");
0067 selectHighPurity_ = iConfig.getParameter<bool>("selectHighPurity");
0068
0069 fillTrackerKink_ = false;
0070 if (fillGlobalTrackQuality_)
0071 fillTrackerKink_ = iConfig.getParameter<bool>("fillTrackerKink");
0072
0073 ptThresholdToFillCandidateP4WithGlobalFit_ =
0074 iConfig.getParameter<double>("ptThresholdToFillCandidateP4WithGlobalFit");
0075 sigmaThresholdToFillCandidateP4WithGlobalFit_ =
0076 iConfig.getParameter<double>("sigmaThresholdToFillCandidateP4WithGlobalFit");
0077 caloCut_ = iConfig.getParameter<double>("minCaloCompatibility");
0078 arbClean_ = iConfig.getParameter<bool>("runArbitrationCleaner");
0079
0080
0081 const edm::ParameterSet parameters = iConfig.getParameter<edm::ParameterSet>("TrackAssociatorParameters");
0082 edm::ConsumesCollector iC = consumesCollector();
0083 parameters_.loadParameters(parameters, iC);
0084
0085
0086 edm::ParameterSet timingParameters = iConfig.getParameter<edm::ParameterSet>("TimingFillerParameters");
0087 theTimingFiller_ = std::make_unique<MuonTimingFiller>(timingParameters, consumesCollector());
0088
0089
0090 if (fillShowerDigis_ && fillMatching_) {
0091 edm::ParameterSet showerDigiParameters = iConfig.getParameter<edm::ParameterSet>("ShowerDigiFillerParameters");
0092 theShowerDigiFiller_ = std::make_unique<MuonShowerDigiFiller>(showerDigiParameters, consumesCollector());
0093 } else {
0094 theShowerDigiFiller_ = std::make_unique<MuonShowerDigiFiller>();
0095 }
0096
0097 if (fillCaloCompatibility_) {
0098
0099 const auto caloParams = iConfig.getParameter<edm::ParameterSet>("MuonCaloCompatibility");
0100 muonCaloCompatibility_.configure(caloParams);
0101 }
0102
0103 if (fillIsolation_) {
0104
0105 edm::ParameterSet caloExtractorPSet = iConfig.getParameter<edm::ParameterSet>("CaloExtractorPSet");
0106 std::string caloExtractorName = caloExtractorPSet.getParameter<std::string>("ComponentName");
0107 muIsoExtractorCalo_ =
0108 IsoDepositExtractorFactory::get()->create(caloExtractorName, caloExtractorPSet, consumesCollector());
0109
0110 edm::ParameterSet trackExtractorPSet = iConfig.getParameter<edm::ParameterSet>("TrackExtractorPSet");
0111 std::string trackExtractorName = trackExtractorPSet.getParameter<std::string>("ComponentName");
0112 muIsoExtractorTrack_ =
0113 IsoDepositExtractorFactory::get()->create(trackExtractorName, trackExtractorPSet, consumesCollector());
0114
0115 edm::ParameterSet jetExtractorPSet = iConfig.getParameter<edm::ParameterSet>("JetExtractorPSet");
0116 std::string jetExtractorName = jetExtractorPSet.getParameter<std::string>("ComponentName");
0117 muIsoExtractorJet_ =
0118 IsoDepositExtractorFactory::get()->create(jetExtractorName, jetExtractorPSet, consumesCollector());
0119 }
0120 if (fillIsolation_ && writeIsoDeposits_) {
0121 trackDepositName_ = iConfig.getParameter<std::string>("trackDepositName");
0122 produces<reco::IsoDepositMap>(trackDepositName_);
0123 ecalDepositName_ = iConfig.getParameter<std::string>("ecalDepositName");
0124 produces<reco::IsoDepositMap>(ecalDepositName_);
0125 hcalDepositName_ = iConfig.getParameter<std::string>("hcalDepositName");
0126 produces<reco::IsoDepositMap>(hcalDepositName_);
0127 hoDepositName_ = iConfig.getParameter<std::string>("hoDepositName");
0128 produces<reco::IsoDepositMap>(hoDepositName_);
0129 jetDepositName_ = iConfig.getParameter<std::string>("jetDepositName");
0130 produces<reco::IsoDepositMap>(jetDepositName_);
0131 }
0132
0133 inputCollectionLabels_ = iConfig.getParameter<std::vector<edm::InputTag> >("inputCollectionLabels");
0134 const auto inputCollectionTypes = iConfig.getParameter<std::vector<std::string> >("inputCollectionTypes");
0135 if (inputCollectionLabels_.size() != inputCollectionTypes.size())
0136 throw cms::Exception("ConfigurationError")
0137 << "Number of input collection labels is different from number of types. "
0138 << "For each collection label there should be exactly one collection type specified.";
0139 if (inputCollectionLabels_.size() > 8 || inputCollectionLabels_.empty())
0140 throw cms::Exception("ConfigurationError") << "Number of input collections should be from 1 to 8.";
0141
0142 debugWithTruthMatching_ = iConfig.getParameter<bool>("debugWithTruthMatching");
0143 if (debugWithTruthMatching_) {
0144 edm::LogWarning("MuonIdentification")
0145 << "========================================================================\n"
0146 << "Debugging mode with truth matching is turned on!!! Make sure you understand what you are doing!\n"
0147 << "========================================================================\n";
0148
0149 globalGeomToken_ = esConsumes();
0150 }
0151 if (fillGlobalTrackQuality_) {
0152 const auto& glbQualTag = iConfig.getParameter<edm::InputTag>("globalTrackQualityInputTag");
0153 glbQualToken_ = consumes<edm::ValueMap<reco::MuonQuality> >(glbQualTag);
0154 }
0155
0156 if (fillTrackerKink_) {
0157 trackerKinkFinder_ =
0158 std::make_unique<MuonKinkFinder>(iConfig.getParameter<edm::ParameterSet>("TrackerKinkFinderParameters"), iC);
0159 }
0160
0161 if (selectHighPurity_) {
0162 const auto& pvTag = iConfig.getParameter<edm::InputTag>("pvInputTag");
0163 pvToken_ = mayConsume<reco::VertexCollection>(pvTag);
0164 }
0165
0166
0167 meshAlgo_ = std::make_unique<MuonMesh>(iConfig.getParameter<edm::ParameterSet>("arbitrationCleanerOptions"));
0168
0169 edm::InputTag rpcHitTag("rpcRecHits");
0170 rpcHitToken_ = consumes<RPCRecHitCollection>(rpcHitTag);
0171
0172 edm::InputTag gemHitTag("gemRecHits");
0173 gemHitToken_ = consumes<GEMRecHitCollection>(gemHitTag);
0174
0175
0176 inputCollectionTypes_.resize(inputCollectionLabels_.size());
0177 for (unsigned int i = 0; i < inputCollectionLabels_.size(); ++i) {
0178 const auto inputLabel = inputCollectionLabels_[i];
0179 const auto inputType = ICTypes::toKey(inputCollectionTypes[i]);
0180
0181 if (inputType == ICTypes::INNER_TRACKS) {
0182 innerTrackCollectionToken_ = consumes<reco::TrackCollection>(inputLabel);
0183 } else if (inputType == ICTypes::OUTER_TRACKS && outerTrackCollectionToken_.isUninitialized()) {
0184 outerTrackCollectionToken_ = consumes<reco::TrackCollection>(inputLabel);
0185 } else if (inputType == ICTypes::OUTER_TRACKS) {
0186 outerTrackSecondaryCollectionToken_ = consumes<reco::TrackCollection>(inputLabel);
0187 } else if (inputType == ICTypes::LINKS) {
0188 linkCollectionToken_ = consumes<reco::MuonTrackLinksCollection>(inputLabel);
0189 } else if (inputType == ICTypes::MUONS) {
0190 muonCollectionToken_ = consumes<reco::MuonCollection>(inputLabel);
0191 } else if (inputType == ICTypes::TEV_FIRSTHIT) {
0192 tpfmsCollectionToken_ = consumes<reco::TrackToTrackMap>(inputLabel);
0193 } else if (fillGlobalTrackRefits_ && inputType == ICTypes::TEV_PICKY) {
0194 pickyCollectionToken_ = consumes<reco::TrackToTrackMap>(inputLabel);
0195 } else if (fillGlobalTrackRefits_ && inputType == ICTypes::TEV_DYT) {
0196 dytCollectionToken_ = consumes<reco::TrackToTrackMap>(inputCollectionLabels_.at(i));
0197 }
0198
0199 inputCollectionTypes_[i] = inputType;
0200 }
0201 }
0202
0203 MuonIdProducer::~MuonIdProducer() {
0204
0205 }
0206
0207 void MuonIdProducer::init(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0208 innerTrackCollectionHandle_.clear();
0209 outerTrackCollectionHandle_.clear();
0210 outerTrackSecondaryCollectionHandle_.clear();
0211 linkCollectionHandle_.clear();
0212 muonCollectionHandle_.clear();
0213
0214 tpfmsCollectionHandle_.clear();
0215 pickyCollectionHandle_.clear();
0216 dytCollectionHandle_.clear();
0217
0218 trackAssociator_.setPropagator(&iSetup.getData(propagatorToken_));
0219
0220 if (fillTrackerKink_)
0221 trackerKinkFinder_->init(iSetup);
0222
0223 for (unsigned int i = 0; i < inputCollectionLabels_.size(); ++i) {
0224 const auto& inputLabel = inputCollectionLabels_[i];
0225 const auto inputType = inputCollectionTypes_[i];
0226 if (inputType == ICTypes::INNER_TRACKS) {
0227 iEvent.getByToken(innerTrackCollectionToken_, innerTrackCollectionHandle_);
0228 if (!innerTrackCollectionHandle_.isValid())
0229 throw cms::Exception("FatalError") << "Failed to get input track collection with label: " << inputLabel;
0230 LogTrace("MuonIdentification") << "Number of input inner tracks: " << innerTrackCollectionHandle_->size();
0231 } else if (inputType == ICTypes::OUTER_TRACKS && !outerTrackCollectionHandle_.isValid()) {
0232 iEvent.getByToken(outerTrackCollectionToken_, outerTrackCollectionHandle_);
0233 if (!outerTrackCollectionHandle_.isValid())
0234 throw cms::Exception("FatalError") << "Failed to get input track collection with label: " << inputLabel;
0235 LogTrace("MuonIdentification") << "Number of input outer tracks: " << outerTrackCollectionHandle_->size();
0236 } else if (inputType == ICTypes::OUTER_TRACKS) {
0237 iEvent.getByToken(outerTrackSecondaryCollectionToken_, outerTrackSecondaryCollectionHandle_);
0238 if (!outerTrackSecondaryCollectionHandle_.isValid())
0239 throw cms::Exception("FatalError") << "Failed to get input track collection with label: " << inputLabel;
0240 LogTrace("MuonIdentification") << "Number of input outer secondary tracks: "
0241 << outerTrackSecondaryCollectionHandle_->size();
0242 } else if (inputType == ICTypes::LINKS) {
0243 iEvent.getByToken(linkCollectionToken_, linkCollectionHandle_);
0244 if (!linkCollectionHandle_.isValid())
0245 throw cms::Exception("FatalError") << "Failed to get input link collection with label: " << inputLabel;
0246 LogTrace("MuonIdentification") << "Number of input links: " << linkCollectionHandle_->size();
0247 } else if (inputType == ICTypes::MUONS) {
0248 iEvent.getByToken(muonCollectionToken_, muonCollectionHandle_);
0249 if (!muonCollectionHandle_.isValid())
0250 throw cms::Exception("FatalError") << "Failed to get input muon collection with label: " << inputLabel;
0251 LogTrace("MuonIdentification") << "Number of input muons: " << muonCollectionHandle_->size();
0252 } else if (fillGlobalTrackRefits_ && inputType == ICTypes::TEV_FIRSTHIT) {
0253 iEvent.getByToken(tpfmsCollectionToken_, tpfmsCollectionHandle_);
0254 if (!tpfmsCollectionHandle_.isValid())
0255 throw cms::Exception("FatalError") << "Failed to get input muon collection with label: " << inputLabel;
0256 LogTrace("MuonIdentification") << "Number of input muons: " << tpfmsCollectionHandle_->size();
0257 } else if (fillGlobalTrackRefits_ && inputType == ICTypes::TEV_PICKY) {
0258 iEvent.getByToken(pickyCollectionToken_, pickyCollectionHandle_);
0259 if (!pickyCollectionHandle_.isValid())
0260 throw cms::Exception("FatalError") << "Failed to get input muon collection with label: " << inputLabel;
0261 LogTrace("MuonIdentification") << "Number of input muons: " << pickyCollectionHandle_->size();
0262 } else if (fillGlobalTrackRefits_ && inputType == ICTypes::TEV_DYT) {
0263 iEvent.getByToken(dytCollectionToken_, dytCollectionHandle_);
0264 if (!dytCollectionHandle_.isValid())
0265 throw cms::Exception("FatalError") << "Failed to get input muon collection with label: " << inputLabel;
0266 LogTrace("MuonIdentification") << "Number of input muons: " << dytCollectionHandle_->size();
0267 } else
0268 throw cms::Exception("FatalError") << "Unknown input collection type: #" << ICTypes::toStr(inputType);
0269 }
0270
0271 iEvent.getByToken(rpcHitToken_, rpcHitHandle_);
0272 iEvent.getByToken(gemHitToken_, gemHitHandle_);
0273 if (fillGlobalTrackQuality_)
0274 iEvent.getByToken(glbQualToken_, glbQualHandle_);
0275 if (selectHighPurity_)
0276 iEvent.getByToken(pvToken_, pvHandle_);
0277 }
0278
0279 reco::Muon MuonIdProducer::makeMuon(edm::Event& iEvent,
0280 const edm::EventSetup& iSetup,
0281 const reco::TrackRef& track,
0282 MuonIdProducer::TrackType type) {
0283 LogTrace("MuonIdentification") << "Creating a muon from a track " << track.get()->pt()
0284 << " Pt (GeV), eta: " << track.get()->eta();
0285 reco::Muon aMuon(makeMuon(*(track.get())));
0286
0287 LogTrace("MuonIdentification") << "Muon created from a track ";
0288
0289 aMuon.setMuonTrack(type, track);
0290 aMuon.setBestTrack(type);
0291 aMuon.setTunePBestTrack(type);
0292
0293 LogTrace("MuonIdentification")
0294 << "Muon created from a track and setMuonBestTrack, setBestTrack and setTunePBestTrack called";
0295
0296 return aMuon;
0297 }
0298
0299 reco::CaloMuon MuonIdProducer::makeCaloMuon(const reco::Muon& muon) {
0300 LogTrace("MuonIdentification") << "Creating a CaloMuon from a Muon";
0301
0302 reco::CaloMuon aMuon;
0303 aMuon.setInnerTrack(muon.innerTrack());
0304
0305 if (muon.isEnergyValid())
0306 aMuon.setCalEnergy(muon.calEnergy());
0307
0308 if (fillCaloCompatibility_)
0309 aMuon.setCaloCompatibility(muonCaloCompatibility_.evaluate(muon));
0310 return aMuon;
0311 }
0312
0313 reco::Muon MuonIdProducer::makeMuon(const reco::MuonTrackLinks& links) {
0314 LogTrace("MuonIdentification") << "Creating a muon from a link to tracks object";
0315
0316 reco::Muon aMuon;
0317 reco::Muon::MuonTrackTypePair chosenTrack;
0318 reco::TrackRef tpfmsRef;
0319 reco::TrackRef pickyRef;
0320 reco::TrackRef dytRef;
0321 bool useSigmaSwitch = false;
0322
0323 if (tpfmsCollectionHandle_.isValid() && !tpfmsCollectionHandle_.failedToGet() && pickyCollectionHandle_.isValid() &&
0324 !pickyCollectionHandle_.failedToGet()) {
0325 tpfmsRef = muon::getTevRefitTrack(links.globalTrack(), *tpfmsCollectionHandle_);
0326 pickyRef = muon::getTevRefitTrack(links.globalTrack(), *pickyCollectionHandle_);
0327 dytRef = muon::getTevRefitTrack(links.globalTrack(), *dytCollectionHandle_);
0328
0329 if (tpfmsRef.isNull() && pickyRef.isNull() && dytRef.isNull()) {
0330 edm::LogWarning("MakeMuonWithTEV") << "Failed to get TEV refits, fall back to sigma switch.";
0331 useSigmaSwitch = true;
0332 }
0333 } else {
0334 useSigmaSwitch = true;
0335 }
0336
0337 if (useSigmaSwitch) {
0338 chosenTrack = muon::sigmaSwitch(links.globalTrack(),
0339 links.trackerTrack(),
0340 sigmaThresholdToFillCandidateP4WithGlobalFit_,
0341 ptThresholdToFillCandidateP4WithGlobalFit_);
0342 } else {
0343 chosenTrack = muon::tevOptimized(links.globalTrack(),
0344 links.trackerTrack(),
0345 tpfmsRef,
0346 pickyRef,
0347 dytRef,
0348 ptThresholdToFillCandidateP4WithGlobalFit_);
0349 }
0350 aMuon = makeMuon(*chosenTrack.first);
0351 aMuon.setInnerTrack(links.trackerTrack());
0352 aMuon.setOuterTrack(links.standAloneTrack());
0353 aMuon.setGlobalTrack(links.globalTrack());
0354 aMuon.setBestTrack(chosenTrack.second);
0355 aMuon.setTunePBestTrack(chosenTrack.second);
0356
0357 if (fillGlobalTrackRefits_) {
0358 if (tpfmsCollectionHandle_.isValid() && !tpfmsCollectionHandle_.failedToGet()) {
0359 reco::TrackToTrackMap::const_iterator it = tpfmsCollectionHandle_->find(links.globalTrack());
0360 if (it != tpfmsCollectionHandle_->end())
0361 aMuon.setMuonTrack(reco::Muon::TPFMS, (it->val));
0362 }
0363 if (pickyCollectionHandle_.isValid() && !pickyCollectionHandle_.failedToGet()) {
0364 reco::TrackToTrackMap::const_iterator it = pickyCollectionHandle_->find(links.globalTrack());
0365 if (it != pickyCollectionHandle_->end())
0366 aMuon.setMuonTrack(reco::Muon::Picky, (it->val));
0367 }
0368 if (dytCollectionHandle_.isValid() && !dytCollectionHandle_.failedToGet()) {
0369 reco::TrackToTrackMap::const_iterator it = dytCollectionHandle_->find(links.globalTrack());
0370 if (it != dytCollectionHandle_->end())
0371 aMuon.setMuonTrack(reco::Muon::DYT, (it->val));
0372 }
0373 }
0374 return aMuon;
0375 }
0376
0377 bool MuonIdProducer::isGoodTrack(const reco::Track& track) {
0378
0379 const double p = track.p();
0380 const double pt = track.pt();
0381 if (pt < minPt_ || (p < minP_ && p < minPCaloMuon_)) {
0382 LogTrace("MuonIdentification") << "Skipped low momentum track (Pt,P): " << pt << ", " << track.p() << " GeV";
0383 return false;
0384 }
0385
0386
0387 const double eta = track.eta();
0388 const double absEta = std::abs(eta);
0389 if (absEta > maxAbsEta_) {
0390 LogTrace("MuonIdentification") << "Skipped track with large pseudo rapidity (Eta: " << track.eta() << " )";
0391 return false;
0392 }
0393
0394 return true;
0395 }
0396
0397 unsigned int MuonIdProducer::chamberId(const DetId& id) {
0398 if (id.det() != DetId::Muon)
0399 return 0;
0400
0401 const auto subdetId = id.subdetId();
0402 if (subdetId == MuonSubdetId::DT) {
0403 return DTChamberId(id.rawId()).rawId();
0404 } else if (subdetId == MuonSubdetId::CSC) {
0405 return CSCDetId(id.rawId()).chamberId().rawId();
0406 }
0407
0408 return 0;
0409 }
0410
0411 int MuonIdProducer::overlap(const reco::Muon& muon, const reco::Track& track) {
0412 if (!muon.isMatchesValid() || track.extra().isNull() || track.extra()->recHitsSize() == 0)
0413 return 0;
0414
0415 int numberOfCommonDetIds = 0;
0416 const std::vector<reco::MuonChamberMatch>& matches(muon.matches());
0417 for (const auto& match : matches) {
0418 if (match.segmentMatches.empty())
0419 continue;
0420
0421 bool foundCommonDetId = false;
0422 for (auto hit = track.extra()->recHitsBegin(); hit != track.extra()->recHitsEnd(); ++hit) {
0423
0424
0425
0426
0427 if (chamberId((*hit)->geographicalId()) == match.id.rawId()) {
0428 foundCommonDetId = true;
0429 break;
0430 }
0431 }
0432 if (foundCommonDetId) {
0433 ++numberOfCommonDetIds;
0434 break;
0435 }
0436 }
0437 return numberOfCommonDetIds;
0438 }
0439
0440 void MuonIdProducer::beginRun(const edm::Run& iRun, const edm::EventSetup& iSetup) {
0441 meshAlgo_->setCSCGeometry(&iSetup.getData(geomTokenRun_));
0442
0443 if (fillShowerDigis_ && fillMatching_)
0444 theShowerDigiFiller_->getES(iSetup);
0445 }
0446
0447 bool validateGlobalMuonPair(const reco::MuonTrackLinks& goodMuon, const reco::MuonTrackLinks& badMuon) {
0448 const int nHitsGood = goodMuon.globalTrack()->hitPattern().numberOfValidMuonHits();
0449 const int nHitsBad = badMuon.globalTrack()->hitPattern().numberOfValidMuonHits();
0450 if (std::min(nHitsGood, nHitsBad) > 10) {
0451 const double chi2Good = goodMuon.globalTrack()->normalizedChi2();
0452 const double chi2Bad = badMuon.globalTrack()->normalizedChi2();
0453 return (chi2Good <= chi2Bad);
0454 }
0455
0456 return (nHitsGood >= nHitsBad);
0457 }
0458
0459 void MuonIdProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0460 auto outputMuons = std::make_unique<reco::MuonCollection>();
0461 auto caloMuons = std::make_unique<reco::CaloMuonCollection>();
0462
0463 init(iEvent, iSetup);
0464
0465 if (fillShowerDigis_ && fillMatching_)
0466 theShowerDigiFiller_->getDigis(iEvent);
0467
0468
0469
0470
0471 if (muonCollectionHandle_.isValid()) {
0472 for (const auto& muon : *muonCollectionHandle_) {
0473 outputMuons->push_back(muon);
0474 }
0475 }
0476
0477
0478 if (linkCollectionHandle_.isValid()) {
0479 const auto nLink = linkCollectionHandle_->size();
0480 std::vector<bool> goodmuons(nLink, true);
0481 if (nLink > 1) {
0482
0483 for (unsigned int i = 0; i < nLink - 1; ++i) {
0484 const auto& iLink = linkCollectionHandle_->at(i);
0485 if (iLink.trackerTrack().isNull() || !checkLinks(&iLink))
0486 continue;
0487 for (unsigned int j = i + 1; j < nLink; ++j) {
0488 const auto& jLink = linkCollectionHandle_->at(j);
0489 if (!checkLinks(&jLink))
0490 continue;
0491 if (iLink.trackerTrack() == jLink.trackerTrack()) {
0492
0493
0494
0495
0496 if (validateGlobalMuonPair(iLink, jLink))
0497 goodmuons[j] = false;
0498 else
0499 goodmuons[i] = false;
0500 }
0501 }
0502 }
0503
0504 for (unsigned int i = 0; i < nLink - 1; ++i) {
0505 if (!goodmuons[i])
0506 continue;
0507 const auto& iLink = linkCollectionHandle_->at(i);
0508 if (iLink.standAloneTrack().isNull() || !checkLinks(&iLink))
0509 continue;
0510 for (unsigned int j = i + 1; j < nLink; ++j) {
0511 if (!goodmuons[j])
0512 continue;
0513 const auto& jLink = linkCollectionHandle_->at(j);
0514 if (!checkLinks(&jLink))
0515 continue;
0516 if (iLink.standAloneTrack() == jLink.standAloneTrack()) {
0517 if (validateGlobalMuonPair(iLink, jLink))
0518 goodmuons[j] = false;
0519 else
0520 goodmuons[i] = false;
0521 }
0522 }
0523 }
0524 }
0525 for (unsigned int i = 0; i < nLink; ++i) {
0526 if (!goodmuons[i])
0527 continue;
0528 const auto& iLink = linkCollectionHandle_->at(i);
0529 if (!checkLinks(&iLink))
0530 continue;
0531
0532 bool newMuon = true;
0533 for (const auto& muon : *outputMuons) {
0534 if (muon.track() == iLink.trackerTrack() && muon.standAloneMuon() == iLink.standAloneTrack() &&
0535 muon.combinedMuon() == iLink.globalTrack()) {
0536 newMuon = false;
0537 break;
0538 }
0539 }
0540 if (newMuon) {
0541 outputMuons->push_back(makeMuon(iLink));
0542 outputMuons->back().setType(reco::Muon::GlobalMuon | reco::Muon::StandAloneMuon);
0543 }
0544 }
0545 }
0546
0547
0548 if (innerTrackCollectionHandle_.isValid()) {
0549 LogTrace("MuonIdentification") << "Creating tracker muons";
0550 std::vector<TrackDetectorAssociator::Direction> directions1, directions2;
0551 directions1.push_back(TrackDetectorAssociator::InsideOut);
0552 directions1.push_back(TrackDetectorAssociator::OutsideIn);
0553 directions2.push_back(TrackDetectorAssociator::Any);
0554
0555 const GlobalTrackingGeometry* geometry = nullptr;
0556 if (debugWithTruthMatching_) {
0557 geometry = &iSetup.getData(globalGeomToken_);
0558 }
0559
0560 for (unsigned int i = 0; i < innerTrackCollectionHandle_->size(); ++i) {
0561 const reco::Track& track = innerTrackCollectionHandle_->at(i);
0562 if (!isGoodTrack(track))
0563 continue;
0564 if (selectHighPurity_ && !track.quality(reco::TrackBase::highPurity)) {
0565 const reco::VertexCollection* recoVertices = pvHandle_.product();
0566 if (!(*recoVertices)[0].isFake())
0567 continue;
0568 }
0569 const auto& trackRef = reco::TrackRef(innerTrackCollectionHandle_, i);
0570 bool splitTrack = false;
0571 if (track.extra().isAvailable() && TrackDetectorAssociator::crossedIP(track))
0572 splitTrack = true;
0573 const auto& directions = splitTrack ? directions1 : directions2;
0574 for (const auto direction : directions) {
0575
0576 reco::Muon trackerMuon(makeMuon(iEvent, iSetup, trackRef, reco::Muon::InnerTrack));
0577 fillMuonId(iEvent, iSetup, trackerMuon, direction);
0578
0579 if (debugWithTruthMatching_) {
0580
0581
0582 MuonIdTruthInfo::truthMatchMuon(iEvent, *geometry, trackerMuon);
0583 }
0584
0585
0586
0587
0588 bool newMuon = true;
0589 const bool goodTrackerMuon = isGoodTrackerMuon(trackerMuon);
0590 const bool goodRPCMuon = isGoodRPCMuon(trackerMuon);
0591 const bool goodGEMMuon = isGoodGEMMuon(trackerMuon);
0592 const bool goodME0Muon = isGoodME0Muon(trackerMuon);
0593 if (goodTrackerMuon)
0594 trackerMuon.setType(trackerMuon.type() | reco::Muon::TrackerMuon);
0595 if (goodRPCMuon)
0596 trackerMuon.setType(trackerMuon.type() | reco::Muon::RPCMuon);
0597 if (goodGEMMuon)
0598 trackerMuon.setType(trackerMuon.type() | reco::Muon::GEMMuon);
0599 if (goodME0Muon)
0600 trackerMuon.setType(trackerMuon.type() | reco::Muon::ME0Muon);
0601
0602 for (auto& muon : *outputMuons) {
0603 if (muon.innerTrack().get() == trackerMuon.innerTrack().get() &&
0604 std::abs(reco::deltaPhi(phiOfMuonInteractionRegion(muon), phiOfMuonInteractionRegion(trackerMuon))) <
0605 M_PI_2) {
0606 newMuon = false;
0607 muon.setMatches(trackerMuon.matches());
0608 if (trackerMuon.isTimeValid())
0609 muon.setTime(trackerMuon.time());
0610 if (trackerMuon.isEnergyValid())
0611 muon.setCalEnergy(trackerMuon.calEnergy());
0612 if (goodTrackerMuon)
0613 muon.setType(muon.type() | reco::Muon::TrackerMuon);
0614 if (goodRPCMuon)
0615 muon.setType(muon.type() | reco::Muon::RPCMuon);
0616 if (goodGEMMuon)
0617 muon.setType(muon.type() | reco::Muon::GEMMuon);
0618 if (goodME0Muon)
0619 muon.setType(muon.type() | reco::Muon::ME0Muon);
0620 LogTrace("MuonIdentification") << "Found a corresponding global muon. Set energy, matches and move on";
0621 break;
0622 }
0623 }
0624 if (newMuon) {
0625 if (goodTrackerMuon || goodRPCMuon || goodGEMMuon || goodME0Muon) {
0626 outputMuons->push_back(trackerMuon);
0627 } else {
0628 LogTrace("MuonIdentification") << "track failed minimal number of muon matches requirement";
0629 const reco::CaloMuon& caloMuon = makeCaloMuon(trackerMuon);
0630 if (isGoodCaloMuon(caloMuon))
0631 caloMuons->push_back(caloMuon);
0632 }
0633 }
0634 }
0635 }
0636 }
0637
0638
0639 if (outerTrackCollectionHandle_.isValid()) {
0640 LogTrace("MuonIdentification") << "Looking for new muons among stand alone muon tracks";
0641 const unsigned int nouter = outerTrackCollectionHandle_->size();
0642 const unsigned int nsecond =
0643 (outerTrackSecondaryCollectionHandle_.isValid()) ? outerTrackSecondaryCollectionHandle_->size() : 0;
0644 for (unsigned int i = 0; i < nouter + nsecond; ++i) {
0645 const auto& outerTrack =
0646 (i < nouter) ? outerTrackCollectionHandle_->at(i) : outerTrackSecondaryCollectionHandle_->at(i - nouter);
0647 reco::TrackRef refToTrack = (i < nouter) ? reco::TrackRef(outerTrackCollectionHandle_, i)
0648 : reco::TrackRef(outerTrackSecondaryCollectionHandle_, i - nouter);
0649
0650
0651 bool newMuon = true;
0652 for (auto& muon : *outputMuons) {
0653 if (!muon.standAloneMuon().isNull()) {
0654
0655 if (muon.standAloneMuon().get() == &outerTrack ||
0656 (muon.standAloneMuon()->extra().isNonnull() &&
0657 muon.standAloneMuon()->extra().get() == outerTrack.extra().get())) {
0658 newMuon = false;
0659 break;
0660 }
0661 } else {
0662
0663
0664
0665
0666
0667
0668
0669 if (overlap(muon, outerTrack) > 0) {
0670 LogTrace("MuonIdentification") << "Found associated tracker muon. Set a reference and move on";
0671 newMuon = false;
0672 muon.setOuterTrack(refToTrack);
0673 muon.setType(muon.type() | reco::Muon::StandAloneMuon);
0674 break;
0675 }
0676 }
0677 }
0678 if (newMuon) {
0679 LogTrace("MuonIdentification") << "No associated stand alone track is found. Making a muon";
0680 outputMuons->push_back(makeMuon(iEvent, iSetup, refToTrack, reco::Muon::OuterTrack));
0681 outputMuons->back().setType(reco::Muon::StandAloneMuon);
0682 }
0683 }
0684 }
0685
0686 if (arbitrateTrackerMuons_) {
0687 fillArbitrationInfo(outputMuons.get(), reco::Muon::TrackerMuon);
0688 fillArbitrationInfo(outputMuons.get(), reco::Muon::GEMMuon);
0689 fillArbitrationInfo(outputMuons.get(), reco::Muon::ME0Muon);
0690 arbitrateMuons(outputMuons.get(), caloMuons.get());
0691 }
0692
0693 LogTrace("MuonIdentification") << "Dress up muons if it's necessary";
0694
0695 const int nMuons = outputMuons->size();
0696
0697 std::vector<reco::MuonTimeExtra> dtTimeColl(nMuons);
0698 std::vector<reco::MuonTimeExtra> cscTimeColl(nMuons);
0699 std::vector<reco::MuonTimeExtra> combinedTimeColl(nMuons);
0700 std::vector<reco::IsoDeposit> trackDepColl(nMuons);
0701 std::vector<reco::IsoDeposit> ecalDepColl(nMuons);
0702 std::vector<reco::IsoDeposit> hcalDepColl(nMuons);
0703 std::vector<reco::IsoDeposit> hoDepColl(nMuons);
0704 std::vector<reco::IsoDeposit> jetDepColl(nMuons);
0705
0706
0707 for (unsigned int i = 0; i < outputMuons->size(); ++i) {
0708 auto& muon = outputMuons->at(i);
0709
0710
0711 if ((fillMatching_ && !muon.isMatchesValid()) || (fillEnergy_ && !muon.isEnergyValid())) {
0712
0713
0714 if (muon.isStandAloneMuon()) {
0715 if (std::abs(reco::deltaPhi(phiOfMuonInteractionRegion(muon), muon.phi())) < M_PI_2) {
0716 fillMuonId(iEvent, iSetup, muon, TrackDetectorAssociator::InsideOut);
0717 } else {
0718 fillMuonId(iEvent, iSetup, muon, TrackDetectorAssociator::OutsideIn);
0719 }
0720 } else {
0721 LogTrace("MuonIdentification") << "THIS SHOULD NEVER HAPPEN";
0722 fillMuonId(iEvent, iSetup, muon);
0723 }
0724 }
0725
0726 if (fillGlobalTrackQuality_) {
0727
0728 fillGlbQuality(iEvent, iSetup, muon);
0729 }
0730 LogDebug("MuonIdentification") << "";
0731
0732 if (fillTrackerKink_) {
0733 fillTrackerKink(muon);
0734 }
0735
0736 if (fillCaloCompatibility_)
0737 muon.setCaloCompatibility(muonCaloCompatibility_.evaluate(muon));
0738
0739 if (fillIsolation_) {
0740 fillMuonIsolation(
0741 iEvent, iSetup, muon, trackDepColl[i], ecalDepColl[i], hcalDepColl[i], hoDepColl[i], jetDepColl[i]);
0742 }
0743
0744
0745 reco::MuonTime muonTime;
0746 reco::MuonTimeExtra dtTime;
0747 reco::MuonTimeExtra cscTime;
0748 reco::MuonTime rpcTime;
0749 reco::MuonTimeExtra combinedTime;
0750
0751 theTimingFiller_->fillTiming(muon, dtTime, cscTime, rpcTime, combinedTime, iEvent, iSetup);
0752
0753 muonTime.nDof = combinedTime.nDof();
0754 muonTime.timeAtIpInOut = combinedTime.timeAtIpInOut();
0755 muonTime.timeAtIpInOutErr = combinedTime.timeAtIpInOutErr();
0756 muonTime.timeAtIpOutIn = combinedTime.timeAtIpOutIn();
0757 muonTime.timeAtIpOutInErr = combinedTime.timeAtIpOutInErr();
0758
0759 muon.setTime(muonTime);
0760 muon.setRPCTime(rpcTime);
0761 dtTimeColl[i] = dtTime;
0762 cscTimeColl[i] = cscTime;
0763 combinedTimeColl[i] = combinedTime;
0764 }
0765
0766 LogTrace("MuonIdentification") << "number of muons produced: " << outputMuons->size();
0767 if (fillMatching_) {
0768 fillArbitrationInfo(outputMuons.get(), reco::Muon::TrackerMuon);
0769 fillArbitrationInfo(outputMuons.get(), reco::Muon::ME0Muon);
0770 fillArbitrationInfo(outputMuons.get(), reco::Muon::GEMMuon);
0771 }
0772 edm::OrphanHandle<reco::MuonCollection> muonHandle = iEvent.put(std::move(outputMuons));
0773
0774 auto fillMap = [](auto refH, auto& vec, edm::Event& ev, const std::string& cAl = "") {
0775 typedef edm::ValueMap<typename std::decay<decltype(vec)>::type::value_type> MapType;
0776 auto oMap = std::make_unique<MapType>();
0777 {
0778 typename MapType::Filler filler(*oMap);
0779 filler.insert(refH, vec.begin(), vec.end());
0780 vec.clear();
0781 filler.fill();
0782 }
0783 ev.put(std::move(oMap), cAl);
0784 };
0785 fillMap(muonHandle, combinedTimeColl, iEvent, "combined");
0786 fillMap(muonHandle, dtTimeColl, iEvent, "dt");
0787 fillMap(muonHandle, cscTimeColl, iEvent, "csc");
0788
0789 if (writeIsoDeposits_ && fillIsolation_) {
0790 fillMap(muonHandle, trackDepColl, iEvent, trackDepositName_);
0791 fillMap(muonHandle, ecalDepColl, iEvent, ecalDepositName_);
0792 fillMap(muonHandle, hcalDepColl, iEvent, hcalDepositName_);
0793 fillMap(muonHandle, hoDepColl, iEvent, hoDepositName_);
0794 fillMap(muonHandle, jetDepColl, iEvent, jetDepositName_);
0795 }
0796
0797 iEvent.put(std::move(caloMuons));
0798 }
0799
0800 bool MuonIdProducer::isGoodTrackerMuon(const reco::Muon& muon) {
0801 if (muon.track()->pt() < minPt_ || muon.track()->p() < minP_)
0802 return false;
0803
0804 int numMatches = muon.numberOfMatches(reco::Muon::NoArbitration);
0805 if (addExtraSoftMuons_ && muon.pt() < 5 && std::abs(muon.eta()) < 1.5 && numMatches >= 1)
0806 return true;
0807 return (numMatches >= minNumberOfMatches_);
0808 }
0809
0810 bool MuonIdProducer::isGoodCaloMuon(const reco::CaloMuon& caloMuon) {
0811 if (!caloMuon.isCaloCompatibilityValid() || caloMuon.caloCompatibility() < caloCut_ || caloMuon.p() < minPCaloMuon_)
0812 return false;
0813 return true;
0814 }
0815
0816 bool MuonIdProducer::isGoodRPCMuon(const reco::Muon& muon) {
0817 if (muon.track()->pt() < minPt_ || muon.track()->p() < minP_)
0818 return false;
0819 if (addExtraSoftMuons_ && muon.pt() < 5 && std::abs(muon.eta()) < 1.5 &&
0820 muon.numberOfMatchedRPCLayers(reco::Muon::RPCHitAndTrackArbitration) > 1)
0821 return true;
0822 return (muon.numberOfMatchedRPCLayers(reco::Muon::RPCHitAndTrackArbitration) > minNumberOfMatches_);
0823 }
0824
0825 bool MuonIdProducer::isGoodGEMMuon(const reco::Muon& muon) {
0826
0827 if (!isGoodTrackerMuon(muon))
0828 return false;
0829 if (muon.track()->pt() < minPt_ || muon.track()->p() < minP_)
0830 return false;
0831
0832 int numMatches = 0;
0833 for (auto& chamberMatch : muon.matches()) {
0834 if (chamberMatch.gemMatches.empty())
0835 continue;
0836 numMatches += chamberMatch.gemMatches.size();
0837 }
0838 return (numMatches + muon.numberOfMatches(reco::Muon::GEMHitAndTrackArbitration)) >= 1;
0839 }
0840
0841 bool MuonIdProducer::isGoodME0Muon(const reco::Muon& muon) {
0842
0843 if (muon.track()->p() < minP_)
0844 return false;
0845 return (muon.numberOfMatches(reco::Muon::ME0SegmentAndTrackArbitration) >= 1);
0846 }
0847
0848 void MuonIdProducer::fillMuonId(edm::Event& iEvent,
0849 const edm::EventSetup& iSetup,
0850 reco::Muon& aMuon,
0851 TrackDetectorAssociator::Direction direction) {
0852 LogTrace("MuonIdentification") << "RecoMuon/MuonIdProducer :: fillMuonId";
0853
0854
0855 const reco::Track* track = nullptr;
0856 if (aMuon.track().isNonnull())
0857 track = aMuon.track().get();
0858 else if (aMuon.standAloneMuon().isNonnull())
0859 track = aMuon.standAloneMuon().get();
0860 else
0861 throw cms::Exception("FatalError")
0862 << "Failed to fill muon id information for a muon with undefined references to tracks";
0863
0864 TrackDetMatchInfo info = trackAssociator_.associate(iEvent, iSetup, *track, parameters_, direction);
0865
0866 LogTrace("MuonIdentification") << "RecoMuon/MuonIdProducer :: fillMuonId :: fillEnergy = " << fillEnergy_;
0867
0868 if (fillEnergy_) {
0869 reco::MuonEnergy muonEnergy;
0870 muonEnergy.em = info.crossedEnergy(TrackDetMatchInfo::EcalRecHits);
0871 muonEnergy.had = info.crossedEnergy(TrackDetMatchInfo::HcalRecHits);
0872 muonEnergy.ho = info.crossedEnergy(TrackDetMatchInfo::HORecHits);
0873 muonEnergy.tower = info.crossedEnergy(TrackDetMatchInfo::TowerTotal);
0874 muonEnergy.emS9 = info.nXnEnergy(TrackDetMatchInfo::EcalRecHits, 1);
0875 muonEnergy.emS25 = info.nXnEnergy(TrackDetMatchInfo::EcalRecHits, 2);
0876 muonEnergy.hadS9 = info.nXnEnergy(TrackDetMatchInfo::HcalRecHits, 1);
0877 muonEnergy.hoS9 = info.nXnEnergy(TrackDetMatchInfo::HORecHits, 1);
0878 muonEnergy.towerS9 = info.nXnEnergy(TrackDetMatchInfo::TowerTotal, 1);
0879 if (storeCrossedHcalRecHits_) {
0880 muonEnergy.crossedHadRecHits.clear();
0881 for (auto hit : info.crossedHcalRecHits) {
0882 reco::HcalMuonRecHit mhit;
0883 mhit.energy = hit->energy();
0884 mhit.chi2 = hit->chi2();
0885 mhit.time = hit->time();
0886 mhit.detId = hit->id();
0887 muonEnergy.crossedHadRecHits.push_back(mhit);
0888 }
0889 }
0890 muonEnergy.ecal_position = info.trkGlobPosAtEcal;
0891 muonEnergy.hcal_position = info.trkGlobPosAtHcal;
0892 if (!info.crossedEcalIds.empty())
0893 muonEnergy.ecal_id = info.crossedEcalIds.front();
0894 if (!info.crossedHcalIds.empty())
0895 muonEnergy.hcal_id = info.crossedHcalIds.front();
0896
0897 DetId emMaxId = info.findMaxDeposition(TrackDetMatchInfo::EcalRecHits, 2);
0898 for (const auto& hit : info.ecalRecHits) {
0899 if (hit->id() != emMaxId)
0900 continue;
0901 muonEnergy.emMax = hit->energy();
0902 muonEnergy.ecal_time = hit->time();
0903 }
0904 DetId hadMaxId = info.findMaxDeposition(TrackDetMatchInfo::HcalRecHits, 1);
0905 for (const auto& hit : info.hcalRecHits) {
0906 if (hit->id() != hadMaxId)
0907 continue;
0908 muonEnergy.hadMax = hit->energy();
0909 muonEnergy.hcal_time = hit->time();
0910 }
0911 aMuon.setCalEnergy(muonEnergy);
0912 }
0913 if (!fillMatching_ && !aMuon.isTrackerMuon() && !aMuon.isRPCMuon() && !aMuon.isGEMMuon())
0914 return;
0915
0916
0917 LogTrace("MuonIdentification") << "RecoMuon/MuonIdProducer :: fillMuonId :: fill muon match info ";
0918 std::vector<reco::MuonChamberMatch> muonChamberMatches;
0919 unsigned int nubmerOfMatchesAccordingToTrackAssociator = 0;
0920 for (const auto& chamber : info.chambers) {
0921 if (chamber.id.subdetId() == MuonSubdetId::RPC && rpcHitHandle_.isValid())
0922 continue;
0923 if (chamber.id.subdetId() == MuonSubdetId::GEM && gemHitHandle_.isValid() &&
0924 GEMDetId(chamber.id.rawId()).station() != 0)
0925 continue;
0926 reco::MuonChamberMatch matchedChamber;
0927
0928 const auto& lErr = chamber.tState.localError();
0929 const auto& lPos = chamber.tState.localPosition();
0930 const auto& lDir = chamber.tState.localDirection();
0931
0932 const auto& localError = lErr.positionError();
0933 matchedChamber.x = lPos.x();
0934 matchedChamber.y = lPos.y();
0935 matchedChamber.xErr = sqrt(localError.xx());
0936 matchedChamber.yErr = sqrt(localError.yy());
0937
0938 matchedChamber.dXdZ = lDir.z() != 0 ? lDir.x() / lDir.z() : 9999;
0939 matchedChamber.dYdZ = lDir.z() != 0 ? lDir.y() / lDir.z() : 9999;
0940
0941 AlgebraicSymMatrix55 trajectoryCovMatrix = lErr.matrix();
0942 matchedChamber.dXdZErr = trajectoryCovMatrix(1, 1) > 0 ? sqrt(trajectoryCovMatrix(1, 1)) : 0;
0943 matchedChamber.dYdZErr = trajectoryCovMatrix(2, 2) > 0 ? sqrt(trajectoryCovMatrix(2, 2)) : 0;
0944
0945 matchedChamber.edgeX = chamber.localDistanceX;
0946 matchedChamber.edgeY = chamber.localDistanceY;
0947
0948 matchedChamber.id = chamber.id;
0949
0950 if (fillShowerDigis_ && fillMatching_) {
0951 theShowerDigiFiller_->fill(matchedChamber);
0952 } else {
0953 theShowerDigiFiller_->fillDefault(matchedChamber);
0954 }
0955
0956 if (!chamber.segments.empty())
0957 ++nubmerOfMatchesAccordingToTrackAssociator;
0958
0959
0960 for (const auto& segment : chamber.segments) {
0961 reco::MuonSegmentMatch matchedSegment;
0962 matchedSegment.x = segment.segmentLocalPosition.x();
0963 matchedSegment.y = segment.segmentLocalPosition.y();
0964 matchedSegment.dXdZ =
0965 segment.segmentLocalDirection.z() ? segment.segmentLocalDirection.x() / segment.segmentLocalDirection.z() : 0;
0966 matchedSegment.dYdZ =
0967 segment.segmentLocalDirection.z() ? segment.segmentLocalDirection.y() / segment.segmentLocalDirection.z() : 0;
0968 matchedSegment.xErr = segment.segmentLocalErrorXX > 0 ? sqrt(segment.segmentLocalErrorXX) : 0;
0969 matchedSegment.yErr = segment.segmentLocalErrorYY > 0 ? sqrt(segment.segmentLocalErrorYY) : 0;
0970 matchedSegment.dXdZErr = segment.segmentLocalErrorDxDz > 0 ? sqrt(segment.segmentLocalErrorDxDz) : 0;
0971 matchedSegment.dYdZErr = segment.segmentLocalErrorDyDz > 0 ? sqrt(segment.segmentLocalErrorDyDz) : 0;
0972 matchedSegment.t0 = segment.t0;
0973 matchedSegment.mask = 0;
0974 matchedSegment.dtSegmentRef = segment.dtSegmentRef;
0975 matchedSegment.cscSegmentRef = segment.cscSegmentRef;
0976 matchedSegment.gemSegmentRef = segment.gemSegmentRef;
0977 matchedSegment.me0SegmentRef = segment.me0SegmentRef;
0978 matchedSegment.hasZed_ = segment.hasZed;
0979 matchedSegment.hasPhi_ = segment.hasPhi;
0980
0981 bool matchedX = false;
0982 bool matchedY = false;
0983 LogTrace("MuonIdentification") << " matching local x, segment x: " << matchedSegment.x
0984 << ", chamber x: " << matchedChamber.x << ", max: " << maxAbsDx_;
0985 LogTrace("MuonIdentification") << " matching local y, segment y: " << matchedSegment.y
0986 << ", chamber y: " << matchedChamber.y << ", max: " << maxAbsDy_;
0987 const double matchedSegChDx = std::abs(matchedSegment.x - matchedChamber.x);
0988 const double matchedSegChDy = std::abs(matchedSegment.y - matchedChamber.y);
0989 if (matchedSegment.xErr > 0 && matchedChamber.xErr > 0)
0990 LogTrace("MuonIdentification") << " xpull: "
0991 << matchedSegChDx / std::sqrt(std::pow(matchedSegment.xErr, 2) +
0992 std::pow(matchedChamber.xErr, 2));
0993 if (matchedSegment.yErr > 0 && matchedChamber.yErr > 0)
0994 LogTrace("MuonIdentification") << " ypull: "
0995 << matchedSegChDy / std::sqrt(std::pow(matchedSegment.yErr, 2) +
0996 std::pow(matchedChamber.yErr, 2));
0997
0998 if (matchedSegChDx < maxAbsDx_)
0999 matchedX = true;
1000 else if (matchedSegment.xErr > 0 && matchedChamber.xErr > 0) {
1001 const double invMatchedSegChPullX2 = std::pow(matchedSegment.xErr, 2) + std::pow(matchedChamber.xErr, 2);
1002 if (matchedSegChDx * matchedSegChDx < maxAbsPullX2_ * invMatchedSegChPullX2)
1003 matchedX = true;
1004 }
1005 if (matchedSegChDy < maxAbsDy_)
1006 matchedY = true;
1007 else if (matchedSegment.yErr > 0 && matchedChamber.yErr > 0) {
1008 const double invMatchedSegChPullY2 = std::pow(matchedSegment.yErr, 2) + std::pow(matchedChamber.yErr, 2);
1009 if (matchedSegChDy * matchedSegChDy < maxAbsPullY2_ * invMatchedSegChPullY2)
1010 matchedY = true;
1011 }
1012 if (matchedX && matchedY) {
1013 if (matchedChamber.id.subdetId() == MuonSubdetId::ME0)
1014 matchedChamber.me0Matches.push_back(matchedSegment);
1015 else if (matchedChamber.id.subdetId() == MuonSubdetId::GEM)
1016 matchedChamber.gemMatches.push_back(matchedSegment);
1017 else
1018 matchedChamber.segmentMatches.push_back(matchedSegment);
1019 }
1020 }
1021 muonChamberMatches.push_back(matchedChamber);
1022 }
1023
1024
1025 LogTrace("MuonIdentification") << "RecoMuon/MuonIdProducer :: fillMuonId :: fill RPC info";
1026 if (rpcHitHandle_.isValid()) {
1027 for (const auto& chamber : info.chambers) {
1028 if (chamber.id.subdetId() != MuonSubdetId::RPC)
1029 continue;
1030 const auto& lErr = chamber.tState.localError();
1031 const auto& lPos = chamber.tState.localPosition();
1032 const auto& lDir = chamber.tState.localDirection();
1033
1034 reco::MuonChamberMatch matchedChamber;
1035
1036 LocalError localError = lErr.positionError();
1037 matchedChamber.x = lPos.x();
1038 matchedChamber.y = lPos.y();
1039 matchedChamber.xErr = sqrt(localError.xx());
1040 matchedChamber.yErr = sqrt(localError.yy());
1041
1042 matchedChamber.dXdZ = lDir.z() != 0 ? lDir.x() / lDir.z() : 9999;
1043 matchedChamber.dYdZ = lDir.z() != 0 ? lDir.y() / lDir.z() : 9999;
1044
1045 AlgebraicSymMatrix55 trajectoryCovMatrix = lErr.matrix();
1046 matchedChamber.dXdZErr = trajectoryCovMatrix(1, 1) > 0 ? sqrt(trajectoryCovMatrix(1, 1)) : 0;
1047 matchedChamber.dYdZErr = trajectoryCovMatrix(2, 2) > 0 ? sqrt(trajectoryCovMatrix(2, 2)) : 0;
1048
1049 matchedChamber.edgeX = chamber.localDistanceX;
1050 matchedChamber.edgeY = chamber.localDistanceY;
1051
1052 theShowerDigiFiller_->fillDefault(matchedChamber);
1053
1054 matchedChamber.id = chamber.id;
1055
1056 for (const auto& rpcRecHit : *rpcHitHandle_) {
1057 reco::MuonRPCHitMatch rpcHitMatch;
1058
1059 if (rpcRecHit.rawId() != chamber.id.rawId())
1060 continue;
1061
1062 rpcHitMatch.x = rpcRecHit.localPosition().x();
1063 rpcHitMatch.mask = 0;
1064 rpcHitMatch.bx = rpcRecHit.BunchX();
1065
1066 const double absDx = std::abs(rpcRecHit.localPosition().x() - chamber.tState.localPosition().x());
1067 if (absDx <= 20 or absDx * absDx <= 16 * localError.xx())
1068 matchedChamber.rpcMatches.push_back(rpcHitMatch);
1069 }
1070
1071 muonChamberMatches.push_back(matchedChamber);
1072 }
1073 }
1074
1075
1076 LogTrace("MuonIdentification") << "RecoMuon/MuonIdProducer :: fillMuonId :: fill GEM info";
1077 if (gemHitHandle_.isValid()) {
1078 for (const auto& chamber : info.chambers) {
1079
1080 if (chamber.id.subdetId() != MuonSubdetId::GEM || GEMDetId(chamber.id.rawId()).station() == 0 ||
1081 GEMDetId(chamber.id.rawId()).layer() == 0)
1082 continue;
1083 const auto& lErr = chamber.tState.localError();
1084 const auto& lPos = chamber.tState.localPosition();
1085 const auto& lDir = chamber.tState.localDirection();
1086
1087 reco::MuonChamberMatch matchedChamber;
1088
1089 LocalError localError = lErr.positionError();
1090 matchedChamber.x = lPos.x();
1091 matchedChamber.y = lPos.y();
1092 matchedChamber.xErr = sqrt(localError.xx());
1093 matchedChamber.yErr = sqrt(localError.yy());
1094
1095 matchedChamber.dXdZ = lDir.z() != 0 ? lDir.x() / lDir.z() : 9999;
1096 matchedChamber.dYdZ = lDir.z() != 0 ? lDir.y() / lDir.z() : 9999;
1097
1098 AlgebraicSymMatrix55 trajectoryCovMatrix = lErr.matrix();
1099 matchedChamber.dXdZErr = trajectoryCovMatrix(1, 1) > 0 ? sqrt(trajectoryCovMatrix(1, 1)) : 0;
1100 matchedChamber.dYdZErr = trajectoryCovMatrix(2, 2) > 0 ? sqrt(trajectoryCovMatrix(2, 2)) : 0;
1101
1102 matchedChamber.edgeX = chamber.localDistanceX;
1103 matchedChamber.edgeY = chamber.localDistanceY;
1104
1105 theShowerDigiFiller_->fillDefault(matchedChamber);
1106
1107 matchedChamber.id = chamber.id;
1108
1109 for (const auto& gemRecHit : *gemHitHandle_) {
1110 reco::MuonGEMHitMatch gemHitMatch;
1111
1112 if (GEMDetId(gemRecHit.gemId().region(),
1113 gemRecHit.gemId().ring(),
1114 gemRecHit.gemId().station(),
1115 gemRecHit.gemId().layer(),
1116 gemRecHit.gemId().chamber(),
1117 0)
1118 .rawId() != chamber.id.rawId())
1119 continue;
1120
1121 gemHitMatch.x = gemRecHit.localPosition().x();
1122 gemHitMatch.mask = 0;
1123 gemHitMatch.bx = gemRecHit.BunchX();
1124
1125 const double absDx = std::abs(gemRecHit.localPosition().x() - chamber.tState.localPosition().x());
1126 if (absDx <= 5 or absDx * absDx <= 16 * localError.xx())
1127 matchedChamber.gemHitMatches.push_back(gemHitMatch);
1128 }
1129
1130 muonChamberMatches.push_back(matchedChamber);
1131 }
1132 }
1133
1134 aMuon.setMatches(muonChamberMatches);
1135
1136 LogTrace("MuonIdentification") << "number of muon chambers: " << aMuon.matches().size() << "\n"
1137 << "number of chambers with segments according to the associator requirements: "
1138 << nubmerOfMatchesAccordingToTrackAssociator;
1139 LogTrace("MuonIdentification") << "number of segment matches with the producer requirements: "
1140 << aMuon.numberOfMatches(reco::Muon::NoArbitration);
1141
1142
1143 }
1144
1145 void MuonIdProducer::arbitrateMuons(reco::MuonCollection* muons, reco::CaloMuonCollection* caloMuons) {
1146 reco::Muon::ArbitrationType arbitration = reco::Muon::SegmentAndTrackArbitration;
1147
1148
1149 for (reco::MuonCollection::iterator muon = muons->begin(); muon != muons->end();) {
1150 if (muon->isTrackerMuon()) {
1151 int numMatches =
1152 muon->numberOfMatches(reco::Muon::GEMSegmentAndTrackArbitration) + muon->numberOfMatches(arbitration);
1153 if (numMatches < minNumberOfMatches_) {
1154
1155
1156
1157
1158 unsigned int mask = reco::Muon::TrackerMuon | reco::Muon::PFMuon | reco::Muon::GEMMuon;
1159 if ((muon->type() & (~mask)) == 0) {
1160 const reco::CaloMuon& caloMuon = makeCaloMuon(*muon);
1161 if (isGoodCaloMuon(caloMuon))
1162 caloMuons->push_back(caloMuon);
1163 muon = muons->erase(muon);
1164 continue;
1165 } else {
1166 muon->setType(muon->type() & (~(reco::Muon::TrackerMuon | reco::Muon::GEMMuon)));
1167 }
1168 }
1169 }
1170 muon++;
1171 }
1172 }
1173
1174 void MuonIdProducer::fillArbitrationInfo(reco::MuonCollection* pOutputMuons, unsigned int muonType) {
1175
1176
1177
1178 std::vector<std::pair<reco::MuonChamberMatch*, reco::MuonSegmentMatch*> > chamberPairs;
1179 std::vector<std::pair<reco::MuonChamberMatch*, reco::MuonSegmentMatch*> > stationPairs;
1180 std::vector<std::pair<reco::MuonChamberMatch*, reco::MuonSegmentMatch*> >
1181 arbitrationPairs;
1182
1183
1184 for (unsigned int muonIndex1 = 0; muonIndex1 < pOutputMuons->size(); ++muonIndex1) {
1185 auto& muon1 = pOutputMuons->at(muonIndex1);
1186
1187 for (auto& chamber1 : muon1.matches()) {
1188
1189 std::vector<reco::MuonSegmentMatch>* segmentMatches1 = getSegmentMatches(chamber1, muonType);
1190
1191 if (segmentMatches1->empty())
1192 continue;
1193 chamberPairs.clear();
1194
1195 for (auto& segment1 : *segmentMatches1) {
1196 chamberPairs.push_back(std::make_pair(&chamber1, &segment1));
1197 if (!segment1.isMask())
1198 {
1199 arbitrationPairs.clear();
1200 arbitrationPairs.push_back(std::make_pair(&chamber1, &segment1));
1201
1202
1203
1204 if (muon1.type() & muonType) {
1205
1206 for (unsigned int muonIndex2 = muonIndex1 + 1; muonIndex2 < pOutputMuons->size(); ++muonIndex2) {
1207 auto& muon2 = pOutputMuons->at(muonIndex2);
1208
1209 if (!(muon2.type() & muonType))
1210 continue;
1211
1212 for (auto& chamber2 : muon2.matches()) {
1213
1214 std::vector<reco::MuonSegmentMatch>* segmentMatches2 = getSegmentMatches(chamber2, muonType);
1215 for (auto& segment2 : *segmentMatches2) {
1216 if (segment2.isMask())
1217 continue;
1218 if (approxEqual(segment2.x, segment1.x) && approxEqual(segment2.y, segment1.y) &&
1219 approxEqual(segment2.dXdZ, segment1.dXdZ) && approxEqual(segment2.dYdZ, segment1.dYdZ) &&
1220 approxEqual(segment2.xErr, segment1.xErr) && approxEqual(segment2.yErr, segment1.yErr) &&
1221 approxEqual(segment2.dXdZErr, segment1.dXdZErr) &&
1222 approxEqual(segment2.dYdZErr, segment1.dYdZErr)) {
1223 arbitrationPairs.push_back(std::make_pair(&chamber2, &segment2));
1224 }
1225 }
1226 }
1227 }
1228 }
1229
1230
1231 if (arbitrationPairs.empty())
1232 continue;
1233 if (arbitrationPairs.size() == 1) {
1234 arbitrationPairs.front().second->setMask(reco::MuonSegmentMatch::BelongsToTrackByDRSlope);
1235 arbitrationPairs.front().second->setMask(reco::MuonSegmentMatch::BelongsToTrackByDXSlope);
1236 arbitrationPairs.front().second->setMask(reco::MuonSegmentMatch::BelongsToTrackByDR);
1237 arbitrationPairs.front().second->setMask(reco::MuonSegmentMatch::BelongsToTrackByDX);
1238 arbitrationPairs.front().second->setMask(reco::MuonSegmentMatch::Arbitrated);
1239 } else {
1240 sort(arbitrationPairs.begin(),
1241 arbitrationPairs.end(),
1242 SortMuonSegmentMatches(reco::MuonSegmentMatch::BelongsToTrackByDRSlope));
1243 arbitrationPairs.front().second->setMask(reco::MuonSegmentMatch::BelongsToTrackByDRSlope);
1244 sort(arbitrationPairs.begin(),
1245 arbitrationPairs.end(),
1246 SortMuonSegmentMatches(reco::MuonSegmentMatch::BelongsToTrackByDXSlope));
1247 arbitrationPairs.front().second->setMask(reco::MuonSegmentMatch::BelongsToTrackByDXSlope);
1248 sort(arbitrationPairs.begin(),
1249 arbitrationPairs.end(),
1250 SortMuonSegmentMatches(reco::MuonSegmentMatch::BelongsToTrackByDR));
1251 arbitrationPairs.front().second->setMask(reco::MuonSegmentMatch::BelongsToTrackByDR);
1252 sort(arbitrationPairs.begin(),
1253 arbitrationPairs.end(),
1254 SortMuonSegmentMatches(reco::MuonSegmentMatch::BelongsToTrackByDX));
1255 arbitrationPairs.front().second->setMask(reco::MuonSegmentMatch::BelongsToTrackByDX);
1256 for (auto& ap : arbitrationPairs) {
1257 ap.second->setMask(reco::MuonSegmentMatch::Arbitrated);
1258 }
1259 }
1260 }
1261
1262
1263 if (muonType == reco::Muon::TrackerMuon && chamber1.id.subdetId() == MuonSubdetId::CSC && arbClean_ &&
1264 CSCDetId(chamber1.id).ring() == 4) {
1265 for (auto& segment2 : chamber1.segmentMatches) {
1266 if (segment1.cscSegmentRef.isNull() || segment2.cscSegmentRef.isNull())
1267 continue;
1268 if (meshAlgo_->isDuplicateOf(segment1.cscSegmentRef, segment2.cscSegmentRef) &&
1269 (segment2.mask & 0x1e0000) && (segment1.mask & 0x1e0000)) {
1270 segment2.setMask(reco::MuonSegmentMatch::BelongsToTrackByME1aClean);
1271
1272 }
1273 }
1274 }
1275
1276 }
1277
1278
1279 if (chamberPairs.empty())
1280 continue;
1281 if (chamberPairs.size() == 1) {
1282 chamberPairs.front().second->setMask(reco::MuonSegmentMatch::BestInChamberByDRSlope);
1283 chamberPairs.front().second->setMask(reco::MuonSegmentMatch::BestInChamberByDXSlope);
1284 chamberPairs.front().second->setMask(reco::MuonSegmentMatch::BestInChamberByDR);
1285 chamberPairs.front().second->setMask(reco::MuonSegmentMatch::BestInChamberByDX);
1286 } else {
1287 sort(chamberPairs.begin(),
1288 chamberPairs.end(),
1289 SortMuonSegmentMatches(reco::MuonSegmentMatch::BestInChamberByDRSlope));
1290 chamberPairs.front().second->setMask(reco::MuonSegmentMatch::BestInChamberByDRSlope);
1291 sort(chamberPairs.begin(),
1292 chamberPairs.end(),
1293 SortMuonSegmentMatches(reco::MuonSegmentMatch::BestInChamberByDXSlope));
1294 chamberPairs.front().second->setMask(reco::MuonSegmentMatch::BestInChamberByDXSlope);
1295 sort(chamberPairs.begin(),
1296 chamberPairs.end(),
1297 SortMuonSegmentMatches(reco::MuonSegmentMatch::BestInChamberByDR));
1298 chamberPairs.front().second->setMask(reco::MuonSegmentMatch::BestInChamberByDR);
1299 sort(chamberPairs.begin(),
1300 chamberPairs.end(),
1301 SortMuonSegmentMatches(reco::MuonSegmentMatch::BestInChamberByDX));
1302 chamberPairs.front().second->setMask(reco::MuonSegmentMatch::BestInChamberByDX);
1303 }
1304 }
1305
1306
1307 for (int stationIndex = 1; stationIndex < 5; ++stationIndex) {
1308 for (int detectorIndex = 1; detectorIndex <= 5;
1309 ++detectorIndex)
1310 {
1311 stationPairs.clear();
1312
1313
1314 for (auto& chamber : muon1.matches()) {
1315 if (!(chamber.station() == stationIndex && chamber.detector() == detectorIndex))
1316 continue;
1317 std::vector<reco::MuonSegmentMatch>* segmentMatches = getSegmentMatches(chamber, muonType);
1318 if (segmentMatches->empty())
1319 continue;
1320
1321 for (auto& segment : *segmentMatches) {
1322 stationPairs.push_back(std::make_pair(&chamber, &segment));
1323 }
1324 }
1325
1326 if (stationPairs.empty())
1327 continue;
1328 if (stationPairs.size() == 1) {
1329 stationPairs.front().second->setMask(reco::MuonSegmentMatch::BestInStationByDRSlope);
1330 stationPairs.front().second->setMask(reco::MuonSegmentMatch::BestInStationByDXSlope);
1331 stationPairs.front().second->setMask(reco::MuonSegmentMatch::BestInStationByDR);
1332 stationPairs.front().second->setMask(reco::MuonSegmentMatch::BestInStationByDX);
1333 } else {
1334 sort(stationPairs.begin(),
1335 stationPairs.end(),
1336 SortMuonSegmentMatches(reco::MuonSegmentMatch::BestInStationByDRSlope));
1337 stationPairs.front().second->setMask(reco::MuonSegmentMatch::BestInStationByDRSlope);
1338 sort(stationPairs.begin(),
1339 stationPairs.end(),
1340 SortMuonSegmentMatches(reco::MuonSegmentMatch::BestInStationByDXSlope));
1341 stationPairs.front().second->setMask(reco::MuonSegmentMatch::BestInStationByDXSlope);
1342 sort(stationPairs.begin(),
1343 stationPairs.end(),
1344 SortMuonSegmentMatches(reco::MuonSegmentMatch::BestInStationByDR));
1345 stationPairs.front().second->setMask(reco::MuonSegmentMatch::BestInStationByDR);
1346 sort(stationPairs.begin(),
1347 stationPairs.end(),
1348 SortMuonSegmentMatches(reco::MuonSegmentMatch::BestInStationByDX));
1349 stationPairs.front().second->setMask(reco::MuonSegmentMatch::BestInStationByDX);
1350 }
1351 }
1352 }
1353
1354 }
1355
1356 if (arbClean_) {
1357
1358 meshAlgo_->clearMesh();
1359 meshAlgo_->runMesh(pOutputMuons);
1360 }
1361 }
1362
1363 void MuonIdProducer::fillMuonIsolation(edm::Event& iEvent,
1364 const edm::EventSetup& iSetup,
1365 reco::Muon& aMuon,
1366 reco::IsoDeposit& trackDep,
1367 reco::IsoDeposit& ecalDep,
1368 reco::IsoDeposit& hcalDep,
1369 reco::IsoDeposit& hoDep,
1370 reco::IsoDeposit& jetDep) {
1371 const reco::Track* track = nullptr;
1372 if (aMuon.track().isNonnull())
1373 track = aMuon.track().get();
1374 else if (aMuon.standAloneMuon().isNonnull())
1375 track = aMuon.standAloneMuon().get();
1376 else
1377 throw cms::Exception("FatalError")
1378 << "Failed to compute muon isolation information for a muon with undefined references to tracks";
1379
1380 reco::MuonIsolation isoR03, isoR05;
1381
1382
1383 reco::IsoDeposit depTrk = muIsoExtractorTrack_->deposit(iEvent, iSetup, *track);
1384 std::vector<reco::IsoDeposit> caloDeps = muIsoExtractorCalo_->deposits(iEvent, iSetup, *track);
1385 reco::IsoDeposit depJet = muIsoExtractorJet_->deposit(iEvent, iSetup, *track);
1386
1387 if (caloDeps.size() != 3) {
1388 LogTrace("MuonIdentification") << "Failed to fill vector of calorimeter isolation deposits!";
1389 return;
1390 }
1391
1392 reco::IsoDeposit depEcal = caloDeps.at(0);
1393 reco::IsoDeposit depHcal = caloDeps.at(1);
1394 reco::IsoDeposit depHo = caloDeps.at(2);
1395
1396
1397 if (writeIsoDeposits_) {
1398 trackDep = depTrk;
1399 ecalDep = depEcal;
1400 hcalDep = depHcal;
1401 hoDep = depHo;
1402 jetDep = depJet;
1403 }
1404
1405 isoR03.sumPt = depTrk.depositWithin(0.3);
1406 isoR03.emEt = depEcal.depositWithin(0.3);
1407 isoR03.hadEt = depHcal.depositWithin(0.3);
1408 isoR03.hoEt = depHo.depositWithin(0.3);
1409 isoR03.nTracks = depTrk.depositAndCountWithin(0.3).second;
1410 isoR03.nJets = depJet.depositAndCountWithin(0.3).second;
1411 isoR03.trackerVetoPt = depTrk.candEnergy();
1412 isoR03.emVetoEt = depEcal.candEnergy();
1413 isoR03.hadVetoEt = depHcal.candEnergy();
1414 isoR03.hoVetoEt = depHo.candEnergy();
1415
1416 isoR05.sumPt = depTrk.depositWithin(0.5);
1417 isoR05.emEt = depEcal.depositWithin(0.5);
1418 isoR05.hadEt = depHcal.depositWithin(0.5);
1419 isoR05.hoEt = depHo.depositWithin(0.5);
1420 isoR05.nTracks = depTrk.depositAndCountWithin(0.5).second;
1421 isoR05.nJets = depJet.depositAndCountWithin(0.5).second;
1422 isoR05.trackerVetoPt = depTrk.candEnergy();
1423 isoR05.emVetoEt = depEcal.candEnergy();
1424 isoR05.hadVetoEt = depHcal.candEnergy();
1425 isoR05.hoVetoEt = depHo.candEnergy();
1426
1427 aMuon.setIsolation(isoR03, isoR05);
1428 }
1429
1430 reco::Muon MuonIdProducer::makeMuon(const reco::Track& track) {
1431 const double energy = std::sqrt(track.p() * track.p() + 0.105658369 * 0.105658369);
1432 const math::XYZTLorentzVector p4(track.px(), track.py(), track.pz(), energy);
1433 return reco::Muon(track.charge(), p4, track.vertex());
1434 }
1435
1436 double MuonIdProducer::sectorPhi(const DetId& id) {
1437 double phi = 0;
1438 if (id.subdetId() == MuonSubdetId::DT) {
1439 DTChamberId muonId(id.rawId());
1440 if (muonId.sector() <= 12)
1441 phi = (muonId.sector() - 1) / 6. * M_PI;
1442 if (muonId.sector() == 13)
1443 phi = 3 / 6. * M_PI;
1444 if (muonId.sector() == 14)
1445 phi = 9 / 6. * M_PI;
1446 }
1447 if (id.subdetId() == MuonSubdetId::CSC) {
1448 CSCDetId muonId(id.rawId());
1449 phi = M_PI / 4 + (muonId.triggerSector() - 1) / 3. * M_PI;
1450 }
1451 if (phi > M_PI)
1452 phi -= 2 * M_PI;
1453 return phi;
1454 }
1455
1456 double MuonIdProducer::phiOfMuonInteractionRegion(const reco::Muon& muon) const {
1457 if (muon.isStandAloneMuon())
1458 return muon.standAloneMuon()->innerPosition().phi();
1459
1460 if (muon.matches().empty()) {
1461 if (muon.innerTrack().isAvailable() && muon.innerTrack()->extra().isAvailable())
1462 return muon.innerTrack()->outerPosition().phi();
1463 else
1464 return muon.phi();
1465 }
1466 return sectorPhi(muon.matches().at(0).id);
1467 }
1468
1469 void MuonIdProducer::fillGlbQuality(edm::Event& iEvent, const edm::EventSetup& iSetup, reco::Muon& aMuon) {
1470 if (aMuon.isGlobalMuon() && glbQualHandle_.isValid() && !glbQualHandle_.failedToGet()) {
1471 aMuon.setCombinedQuality((*glbQualHandle_)[aMuon.combinedMuon()]);
1472 }
1473
1474 LogDebug("MuonIdentification") << "tkChiVal " << aMuon.combinedQuality().trkRelChi2;
1475 }
1476
1477 void MuonIdProducer::fillTrackerKink(reco::Muon& aMuon) {
1478
1479 if (aMuon.innerTrack().isNull())
1480 return;
1481
1482 reco::MuonQuality quality = (aMuon.isQualityValid() ? aMuon.combinedQuality() : reco::MuonQuality());
1483
1484 const bool filled = trackerKinkFinder_->fillTrkKink(quality, *aMuon.innerTrack());
1485
1486 if (filled || aMuon.isQualityValid())
1487 aMuon.setCombinedQuality(quality);
1488 }
1489
1490 bool MuonIdProducer::checkLinks(const reco::MuonTrackLinks* links) const {
1491 const bool trackBAD = links->trackerTrack().isNull();
1492 const bool staBAD = links->standAloneTrack().isNull();
1493 const bool glbBAD = links->globalTrack().isNull();
1494 if (trackBAD || staBAD || glbBAD) {
1495 edm::LogWarning("muonIDbadLinks") << "Global muon links to constituent tracks are invalid: trkBad " << trackBAD
1496 << " standaloneBad " << staBAD << " globalBad " << glbBAD
1497 << ". There should be no such object. Muon is skipped.";
1498 return false;
1499 }
1500 return true;
1501 }
1502
1503 void MuonIdProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
1504 edm::ParameterSetDescription desc;
1505 desc.setAllowAnything();
1506
1507 desc.add<bool>("arbitrateTrackerMuons", false);
1508 desc.add<bool>("storeCrossedHcalRecHits", false);
1509 desc.add<bool>("fillShowerDigis", false);
1510 desc.ifValue(
1511 edm::ParameterDescription<bool>("selectHighPurity", false, true),
1512 true >> (edm::ParameterDescription<edm::InputTag>("pvInputTag", edm::InputTag("offlinePrimaryVertices"), true)) or
1513 false >> (edm::ParameterDescription<edm::InputTag>("pvInputTag", edm::InputTag(""), true)));
1514
1515 edm::ParameterSetDescription descTrkAsoPar;
1516 descTrkAsoPar.add<edm::InputTag>("GEMSegmentCollectionLabel", edm::InputTag("gemSegments"));
1517 descTrkAsoPar.add<edm::InputTag>("ME0SegmentCollectionLabel", edm::InputTag("me0Segments"));
1518 descTrkAsoPar.add<bool>("useGEM", false);
1519 descTrkAsoPar.add<bool>("useME0", false);
1520 descTrkAsoPar.add<bool>("preselectMuonTracks", false);
1521 descTrkAsoPar.add<edm::InputTag>("RPCHitCollectionLabel", edm::InputTag("rpcRecHits"));
1522 descTrkAsoPar.add<edm::InputTag>("GEMHitCollectionLabel", edm::InputTag("gemRecHits"));
1523 descTrkAsoPar.add<edm::InputTag>("ME0HitCollectionLabel", edm::InputTag("me0RecHits"));
1524 descTrkAsoPar.setAllowAnything();
1525 desc.add<edm::ParameterSetDescription>("TrackAssociatorParameters", descTrkAsoPar);
1526
1527 edm::ParameterSetDescription descJet;
1528 descJet.setAllowAnything();
1529 descJet.add<edm::ParameterSetDescription>("TrackAssociatorParameters", descTrkAsoPar);
1530 desc.add<edm::ParameterSetDescription>("JetExtractorPSet", descJet);
1531
1532 edm::ParameterSetDescription descCalo;
1533 descCalo.setAllowAnything();
1534 descCalo.add<edm::ParameterSetDescription>("TrackAssociatorParameters", descTrkAsoPar);
1535 descCalo.add<bool>("UseEcalRecHitsFlag", false);
1536 descCalo.add<bool>("UseHcalRecHitsFlag", false);
1537 descCalo.add<bool>("UseHORecHitsFlag", false);
1538 descCalo.add<bool>("EcalRecHitThresh", false);
1539 descCalo.add<bool>("HcalCutsFromDB", false);
1540 descCalo.add<int>("MaxSeverityHB", 9);
1541 descCalo.add<int>("MaxSeverityHE", 9);
1542 desc.add<edm::ParameterSetDescription>("CaloExtractorPSet", descCalo);
1543
1544 descriptions.addDefault(desc);
1545 }