Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-01-27 00:42:22 

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