Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:27:01 

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() > 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   //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 && 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   // TimingReport::current()->dump(std::cout);
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   // get calo compatibility
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   // Pt and absolute momentum requirement
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   // Eta requirement
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       // LogTrace("MuonIdentification") << "hit DetId: " << std::hex << hit->get()->geographicalId().rawId() <<
0424       //  "\t hit chamber DetId: " << getChamberId(hit->get()->geographicalId()) <<
0425       //  "\t segment DetId: " << match->id.rawId() << std::dec;
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   // loop over input collections
0469 
0470   // muons first - no cleaning, take as is.
0471   if (muonCollectionHandle_.isValid()) {
0472     for (const auto& muon : *muonCollectionHandle_) {
0473       outputMuons->push_back(muon);
0474     }
0475   }
0476 
0477   // links second ( assume global muon type )
0478   if (linkCollectionHandle_.isValid()) {
0479     const auto nLink = linkCollectionHandle_->size();
0480     std::vector<bool> goodmuons(nLink, true);
0481     if (nLink > 1) {
0482       // check for shared tracker tracks
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             // Tracker track is the essential part that dominates muon resolution
0493             // so taking either muon is fine. All that is important is to preserve
0494             // the muon identification information. If number of hits is small,
0495             // keep the one with large number of hits, otherwise take the smalest chi2/ndof
0496             if (validateGlobalMuonPair(iLink, jLink))
0497               goodmuons[j] = false;
0498             else
0499               goodmuons[i] = false;
0500           }
0501         }
0502       }
0503       // check for shared stand-alone muons.
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       // check if this muon is already in the list
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   // tracker and calo muons are next
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         // make muon
0576         reco::Muon trackerMuon(makeMuon(iEvent, iSetup, trackRef, reco::Muon::InnerTrack));
0577         fillMuonId(iEvent, iSetup, trackerMuon, direction);
0578 
0579         if (debugWithTruthMatching_) {
0580           // add MC hits to a list of matched segments.
0581           // Since it's debugging mode - code is slow
0582           MuonIdTruthInfo::truthMatchMuon(iEvent, *geometry, trackerMuon);
0583         }
0584 
0585         // check if this muon is already in the list
0586         // have to check where muon hits are really located
0587         // to match properly
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   // and at last the stand alone muons
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       // check if this muon is already in the list of global muons
0651       bool newMuon = true;
0652       for (auto& muon : *outputMuons) {
0653         if (!muon.standAloneMuon().isNull()) {
0654           // global muon
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           // tracker muon - no direct links to the standalone muon
0663           // since we have only a few real muons in an event, matching
0664           // the stand alone muon to the tracker muon by DetIds should
0665           // be good enough for association. At the end it's up to a
0666           // user to redefine the association and what it means. Here
0667           // we would like to avoid obvious double counting and we
0668           // tolerate a potential miss association
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   // Fill various information
0707   for (unsigned int i = 0; i < outputMuons->size(); ++i) {
0708     auto& muon = outputMuons->at(i);
0709 
0710     // Fill muonID
0711     if ((fillMatching_ && !muon.isMatchesValid()) || (fillEnergy_ && !muon.isEnergyValid())) {
0712       // predict direction based on the muon interaction region location
0713       // if it's available
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       // Fill global quality information
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     // fill timing information
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   // NoArbitration checks for CSC/DT segments only, also use ME0 segments
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   // require GEMMuon to be a TrackerMuon
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   // need to update min cuts on pt
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   // perform track - detector association
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);    // 3x3 energy
0875     muonEnergy.emS25 = info.nXnEnergy(TrackDetMatchInfo::EcalRecHits, 2);   // 5x5 energy
0876     muonEnergy.hadS9 = info.nXnEnergy(TrackDetMatchInfo::HcalRecHits, 1);   // 3x3 energy
0877     muonEnergy.hoS9 = info.nXnEnergy(TrackDetMatchInfo::HORecHits, 1);      // 3x3 energy
0878     muonEnergy.towerS9 = info.nXnEnergy(TrackDetMatchInfo::TowerTotal, 1);  // 3x3 energy
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     // find maximal energy depositions and their time
0897     DetId emMaxId = info.findMaxDeposition(TrackDetMatchInfo::EcalRecHits, 2);  // max energy deposit in 5x5 shape
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);  // max energy deposit in 3x3 shape
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   // fill muon match info
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;  // Skip RPC chambers, they are taken care of below)
0923     if (chamber.id.subdetId() == MuonSubdetId::GEM && gemHitHandle_.isValid() &&
0924         GEMDetId(chamber.id.rawId()).station() != 0)
0925       continue;  // Skip GE1/1 and 2/1 chambers, they are taken care of below)
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     // DANGEROUS - compiler cannot guaranty parameters ordering
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     // fill segments
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       // test segment
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   // Fill RPC info
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;  // Consider RPC chambers only
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       // DANGEROUS - compiler cannot guaranty parameters ordering
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   // Fill GEM info
1076   LogTrace("MuonIdentification") << "RecoMuon/MuonIdProducer :: fillMuonId :: fill GEM info";
1077   if (gemHitHandle_.isValid()) {
1078     for (const auto& chamber : info.chambers) {
1079       // only GE1/1 and 2/1 are for rechits, reject station 0 and segments (layer==0 for GEMSegment)
1080       if (chamber.id.subdetId() != MuonSubdetId::GEM || GEMDetId(chamber.id.rawId()).station() == 0 ||
1081           GEMDetId(chamber.id.rawId()).layer() == 0)
1082         continue;  // Consider GEM chambers only
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       // DANGEROUS - compiler cannot guaranty parameters ordering
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   // fillTime( iEvent, iSetup, aMuon );
1143 }
1144 
1145 void MuonIdProducer::arbitrateMuons(reco::MuonCollection* muons, reco::CaloMuonCollection* caloMuons) {
1146   reco::Muon::ArbitrationType arbitration = reco::Muon::SegmentAndTrackArbitration;
1147   // arbitrate TrackerMuons
1148   // if a muon was exclusively TrackerMuon check if it can be a calo muon
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         // TrackerMuon failed arbitration
1155         // If not any other base type - erase the element
1156         // (PFMuon is not a base type)
1157         // GEMMuon should be a subset of TrackerMuon, so don't count it either
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   // apply segment flags
1177   //
1178   std::vector<std::pair<reco::MuonChamberMatch*, reco::MuonSegmentMatch*> > chamberPairs;  // for chamber segment sorting
1179   std::vector<std::pair<reco::MuonChamberMatch*, reco::MuonSegmentMatch*> > stationPairs;  // for station segment sorting
1180   std::vector<std::pair<reco::MuonChamberMatch*, reco::MuonSegmentMatch*> >
1181       arbitrationPairs;  // for muon segment arbitration
1182 
1183   // muonIndex1
1184   for (unsigned int muonIndex1 = 0; muonIndex1 < pOutputMuons->size(); ++muonIndex1) {
1185     auto& muon1 = pOutputMuons->at(muonIndex1);
1186     // chamberIter1
1187     for (auto& chamber1 : muon1.matches()) {
1188       // segmentIter1
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())  // has not yet been arbitrated
1198         {
1199           arbitrationPairs.clear();
1200           arbitrationPairs.push_back(std::make_pair(&chamber1, &segment1));
1201 
1202           // find identical segments with which to arbitrate
1203           // tracker muons only
1204           if (muon1.type() & muonType) {
1205             // muonIndex2
1206             for (unsigned int muonIndex2 = muonIndex1 + 1; muonIndex2 < pOutputMuons->size(); ++muonIndex2) {
1207               auto& muon2 = pOutputMuons->at(muonIndex2);
1208               // tracker muons only
1209               if (!(muon2.type() & muonType))
1210                 continue;
1211               // chamberIter2
1212               for (auto& chamber2 : muon2.matches()) {
1213                 // segmentIter2
1214                 std::vector<reco::MuonSegmentMatch>* segmentMatches2 = getSegmentMatches(chamber2, muonType);
1215                 for (auto& segment2 : *segmentMatches2) {
1216                   if (segment2.isMask())
1217                     continue;  // has already been arbitrated
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                 }  // segmentIter2
1226               }    // chamberIter2
1227             }      // muonIndex2
1228           }
1229 
1230           // arbitration segment sort
1231           if (arbitrationPairs.empty())
1232             continue;  // this should never happen
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         // setup me1a cleaning for later
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               //if the track has lost the segment already through normal arbitration no need to do it again.
1272             }
1273           }
1274         }  // mark all ME1/a duplicates that this track owns
1275 
1276       }  // segmentIter1
1277 
1278       // chamber segment sort
1279       if (chamberPairs.empty())
1280         continue;  // this should never happen
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     }  // chamberIter1
1305 
1306     // station segment sort
1307     for (int stationIndex = 1; stationIndex < 5; ++stationIndex) {
1308       for (int detectorIndex = 1; detectorIndex <= 5;
1309            ++detectorIndex)  // 1-5, as in DataFormats/MuonDetId/interface/MuonSubdetId.h
1310       {
1311         stationPairs.clear();
1312 
1313         // chamberIter
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         }  // chamberIter
1325 
1326         if (stationPairs.empty())
1327           continue;  // this may very well happen
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   }  // muonIndex1
1355 
1356   if (arbClean_) {
1357     // clear old mesh, create and prune new mesh!
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   // get deposits
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   //no need to copy outside if we don't write them
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) {  // 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) {  // 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   // the rest is tracker muon only
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();  // makes little sense, but what else can I use
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   // skip muons with no tracks
1479   if (aMuon.innerTrack().isNull())
1480     return;
1481   // get quality from muon if already there, otherwise make empty one
1482   reco::MuonQuality quality = (aMuon.isQualityValid() ? aMuon.combinedQuality() : reco::MuonQuality());
1483   // fill it
1484   const bool filled = trackerKinkFinder_->fillTrkKink(quality, *aMuon.innerTrack());
1485   // if quality was there, or if we filled it, commit to the muon
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   desc.add<edm::ParameterSetDescription>("CaloExtractorPSet", descCalo);
1536 
1537   descriptions.addDefault(desc);
1538 }
[ 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!