Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-09-10 01:03:14 

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

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!