Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-05-15 04:21:52

0001 /*
0002  * CandidateSimMuonMatcher.cc
0003  *
0004  *  Created on: Dec 14, 2020
0005  *      Author: kbunkow
0006  */
0007 
0008 #include "L1Trigger/L1TMuonOverlapPhase1/interface/Tools/CandidateSimMuonMatcher.h"
0009 #include "L1Trigger/L1TMuon/interface/MicroGMTConfiguration.h"
0010 
0011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0012 #include "FWCore/Framework/interface/Event.h"
0013 #include "FWCore/Framework/interface/EventSetup.h"
0014 
0015 #include "DataFormats/GeometrySurface/interface/BoundCylinder.h"
0016 #include "DataFormats/GeometrySurface/interface/SimpleCylinderBounds.h"
0017 #include "DataFormats/GeometrySurface/interface/BoundDisk.h"
0018 #include "DataFormats/GeometrySurface/interface/SimpleDiskBounds.h"
0019 #include "DataFormats/MuonDetId/interface/RPCDetId.h"
0020 
0021 #include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h"
0022 
0023 #include "boost/dynamic_bitset.hpp"
0024 
0025 #include "TFile.h"
0026 #include "TH1D.h"
0027 
0028 double hwGmtPhiToGlobalPhi(int phi) {
0029   double phiGmtUnit = 2. * M_PI / 576.;
0030   return phi * phiGmtUnit;
0031 }
0032 
0033 double foldPhi(double phi) {
0034   if (phi > M_PI)
0035     return (phi - 2 * M_PI);
0036   else if (phi < -M_PI)
0037     return (phi + 2 * M_PI);
0038 
0039   return phi;
0040 }
0041 
0042 CandidateSimMuonMatcher::CandidateSimMuonMatcher(
0043     const edm::ParameterSet& edmCfg,
0044     const OMTFConfiguration* omtfConfig,
0045     const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord>& magneticFieldEsToken,
0046     const edm::ESGetToken<Propagator, TrackingComponentsRecord>& propagatorEsToken)
0047     : omtfConfig(omtfConfig),
0048       edmCfg(edmCfg),
0049       magneticFieldEsToken(magneticFieldEsToken),
0050       propagatorEsToken(propagatorEsToken) {
0051   std::string muonMatcherFileName = edmCfg.getParameter<edm::FileInPath>("muonMatcherFile").fullPath();
0052   TFile inFile(muonMatcherFileName.c_str());
0053   edm::LogImportant("l1tOmtfEventPrint") << " CandidateSimMuonMatcher: using muonMatcherFileName "
0054                                          << muonMatcherFileName << std::endl;
0055 
0056   if (edmCfg.exists("candidateSimMuonMatcherType")) {
0057     if (edmCfg.getParameter<std::string>("candidateSimMuonMatcherType") == "propagation")
0058       usePropagation = true;
0059     else if (edmCfg.getParameter<std::string>("candidateSimMuonMatcherType") == "matchSimple")
0060       usePropagation = false;
0061 
0062     edm::LogImportant("l1tOmtfEventPrint")
0063         << " CandidateSimMuonMatcher: candidateSimMuonMatcherType "
0064         << edmCfg.getParameter<std::string>("candidateSimMuonMatcherType") << std::endl;
0065   }
0066 
0067   edm::LogImportant("l1tOmtfEventPrint") << " CandidateSimMuonMatcher: usePropagation " << usePropagation << std::endl;
0068 
0069   deltaPhiPropCandMean = (TH1D*)inFile.Get("deltaPhiPropCandMean");
0070   deltaPhiPropCandStdDev = (TH1D*)inFile.Get("deltaPhiPropCandStdDev");
0071 }
0072 
0073 CandidateSimMuonMatcher::~CandidateSimMuonMatcher() {}
0074 
0075 void CandidateSimMuonMatcher::beginRun(const edm::EventSetup& eventSetup) {
0076   //TODO use edm::ESWatcher<MagneticField> magneticFieldRecordWatcher;
0077   magField = eventSetup.getHandle(magneticFieldEsToken);
0078   propagator = eventSetup.getHandle(propagatorEsToken);
0079 }
0080 
0081 void CandidateSimMuonMatcher::observeEventBegin(const edm::Event& event) { gbCandidates.clear(); }
0082 
0083 void CandidateSimMuonMatcher::observeProcesorEmulation(unsigned int iProcessor,
0084                                                        l1t::tftype mtfType,
0085                                                        const std::shared_ptr<OMTFinput>& input,
0086                                                        const AlgoMuons& algoCandidates,
0087                                                        const AlgoMuons& gbCandidates,
0088                                                        const std::vector<l1t::RegionalMuonCand>& candMuons) {
0089   //debug
0090   unsigned int procIndx = omtfConfig->getProcIndx(iProcessor, mtfType);
0091   for (auto& gbCandidate : gbCandidates) {
0092     if (gbCandidate->getPtConstr() > 0) {
0093       LogTrace("l1tOmtfEventPrint") << "CandidateSimMuonMatcher::observeProcesorEmulation procIndx" << procIndx << " "
0094                                     << *gbCandidate << std::endl;
0095       this->gbCandidates.emplace_back(gbCandidate);
0096     }
0097   }
0098 }
0099 
0100 bool simTrackIsMuonInOmtf(const SimTrack& simTrack) {
0101   if (std::abs(simTrack.type()) == 13 ||
0102       std::abs(simTrack.type()) == 1000015) {  // 1000015 is stau, todo use other selection (e.g. pt>20) if needed
0103     //only muons
0104   } else
0105     return false;
0106 
0107   //in the overlap, the propagation of muons with pt less then ~3.2 fails - the actual threshold depends slightly on eta,
0108   if (simTrack.momentum().pt() < 2.5)
0109     return false;
0110 
0111   LogTrace("l1tOmtfEventPrint") << "simTrackIsMuonInOmtf: simTrack type " << std::setw(3) << simTrack.type() << " pt "
0112                                 << std::setw(9) << simTrack.momentum().pt() << " eta " << std::setw(9)
0113                                 << simTrack.momentum().eta() << " phi " << std::setw(9) << simTrack.momentum().phi()
0114                                 << std::endl;
0115 
0116   //some margin for matching must be used on top of actual OMTF region,
0117   //i.e. (0.82-1.24)=>(0.72-1.3),
0118   //otherwise many candidates are marked as ghosts
0119   if ((std::abs(simTrack.momentum().eta()) >= 0.72) && (std::abs(simTrack.momentum().eta()) <= 1.3)) {
0120     LogTrace("l1tOmtfEventPrint") << "simTrackIsMuonInOmtf: is in OMTF";
0121   } else {
0122     LogTrace("l1tOmtfEventPrint") << "simTrackIsMuonInOmtf: not in OMTF";
0123     return false;
0124   }
0125 
0126   return true;
0127 }
0128 
0129 bool simTrackIsMuonInOmtfBx0(const SimTrack& simTrack) {
0130   if (simTrack.eventId().bunchCrossing() != 0)
0131     return false;
0132 
0133   return simTrackIsMuonInOmtf(simTrack);
0134 }
0135 
0136 bool simTrackIsMuonInBx0(const SimTrack& simTrack) {
0137   if (std::abs(simTrack.type()) == 13 ||
0138       std::abs(simTrack.type()) == 1000015) {  // 1000015 is stau, todo use other selection (e.g. pt>20) if needed
0139     //only muons
0140     if (simTrack.eventId().bunchCrossing() == 0)
0141       return true;
0142   }
0143   return false;
0144 }
0145 
0146 bool trackingParticleIsMuonInBx0(const TrackingParticle& trackingParticle) {
0147   if (std::abs(trackingParticle.pdgId()) == 13 ||
0148       std::abs(trackingParticle.pdgId()) ==
0149           1000015) {  // 1000015 is stau, todo use other selection (e.g. pt>20) if needed
0150     //only muons
0151     if (trackingParticle.eventId().bunchCrossing() == 0)
0152       return true;
0153   }
0154   return false;
0155 }
0156 
0157 bool trackingParticleIsMuonInOmtfBx0(const TrackingParticle& trackingParticle) {
0158   if (trackingParticle.eventId().bunchCrossing() != 0)
0159     return false;
0160 
0161   if (std::abs(trackingParticle.pdgId()) == 13 || std::abs(trackingParticle.pdgId()) == 1000015) {
0162     // 1000015 is stau, todo use other selection (e.g. pt>20) if needed
0163     //only muons
0164   } else
0165     return false;
0166 
0167   //in the overlap, the propagation of muons with pt less then ~3.2 fails, the actual threshold depends slightly on eta,
0168   if (trackingParticle.pt() < 2.5)
0169     return false;
0170 
0171   if (trackingParticle.parentVertex().isNonnull())
0172     LogTrace("l1tOmtfEventPrint") << "trackingParticleIsMuonInOmtfBx0, pdgId " << std::setw(3)
0173                                   << trackingParticle.pdgId() << " pt " << std::setw(9) << trackingParticle.pt()
0174                                   << " eta " << std::setw(9) << trackingParticle.momentum().eta() << " phi "
0175                                   << std::setw(9) << trackingParticle.momentum().phi() << " event "
0176                                   << trackingParticle.eventId().event() << " trackId "
0177                                   << trackingParticle.g4Tracks().at(0).trackId() << " parentVertex Rho "
0178                                   << trackingParticle.parentVertex()->position().Rho() << " eta "
0179                                   << trackingParticle.parentVertex()->position().eta() << " phi "
0180                                   << trackingParticle.parentVertex()->position().phi() << std::endl;
0181   else
0182     LogTrace("l1tOmtfEventPrint") << "trackingParticleIsMuonInOmtfBx0, pdgId " << std::setw(3)
0183                                   << trackingParticle.pdgId() << " pt " << std::setw(9) << trackingParticle.pt()
0184                                   << " eta " << std::setw(9) << trackingParticle.momentum().eta() << " phi "
0185                                   << std::setw(9) << trackingParticle.momentum().phi() << " trackId "
0186                                   << trackingParticle.g4Tracks().at(0).trackId();
0187 
0188   //some margin for matching must be used on top of actual OMTF region,
0189   //i.e. (0.82-1.24)=>(0.72-1.3),
0190   //otherwise many candidates are marked as ghosts
0191   if ((std::abs(trackingParticle.momentum().eta()) >= 0.72) && (std::abs(trackingParticle.momentum().eta()) <= 1.3)) {
0192   } else
0193     return false;
0194 
0195   return true;
0196 }
0197 
0198 bool trackingParticleIsMuonInOmtfEvent0(const TrackingParticle& trackingParticle) {
0199   if (trackingParticle.eventId().event() != 0)
0200     return false;
0201 
0202   return trackingParticleIsMuonInOmtfBx0(trackingParticle);
0203 }
0204 
0205 void CandidateSimMuonMatcher::observeEventEnd(const edm::Event& event,
0206                                               std::unique_ptr<l1t::RegionalMuonCandBxCollection>& finalCandidates) {
0207   LogTrace("l1tOmtfEventPrint") << "\nCandidateSimMuonMatcher::observeEventEnd" << std::endl;
0208   AlgoMuons ghostBustedProcMuons;
0209   std::vector<const l1t::RegionalMuonCand*> ghostBustedRegionalCands =
0210       CandidateSimMuonMatcher::ghostBust(finalCandidates.get(), gbCandidates, ghostBustedProcMuons);
0211 
0212   matchingResults.clear();
0213   if (edmCfg.exists("simTracksTag")) {
0214     edm::Handle<edm::SimTrackContainer> simTraksHandle;
0215     event.getByLabel(edmCfg.getParameter<edm::InputTag>("simTracksTag"), simTraksHandle);
0216 
0217     edm::Handle<edm::SimVertexContainer> simVertices;
0218     event.getByLabel(edmCfg.getParameter<edm::InputTag>("simVertexesTag"), simVertices);
0219 
0220     LogTrace("l1tOmtfEventPrint") << "simTraksHandle size " << simTraksHandle.product()->size() << std::endl;
0221     LogTrace("l1tOmtfEventPrint") << "simVertices size " << simVertices.product()->size() << std::endl;
0222 
0223     if (usePropagation) {
0224       //TODO  use other simTrackFilter if needed  <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
0225       //we dont want to check the eta of the generated muon, as it is on the vertex,
0226       //instead inside match, we check the eta of the propagated track to the second muons station
0227       std::function<bool(const SimTrack&)> const& simTrackFilter = simTrackIsMuonInBx0;  //simTrackIsMuonInOmtfBx0;
0228 
0229       matchingResults = match(ghostBustedRegionalCands,
0230                               ghostBustedProcMuons,
0231                               simTraksHandle.product(),
0232                               simVertices.product(),
0233                               simTrackFilter);
0234     } else {
0235       std::function<bool(const SimTrack&)> const& simTrackFilter = simTrackIsMuonInOmtfBx0;
0236       //simTrackIsMuonInOmtfBx0 provides appropriate eta cut
0237 
0238       matchingResults = matchSimple(ghostBustedRegionalCands,
0239                                     ghostBustedProcMuons,
0240                                     simTraksHandle.product(),
0241                                     simVertices.product(),
0242                                     simTrackFilter);
0243     }
0244 
0245   } else if (edmCfg.exists("trackingParticleTag")) {
0246     edm::Handle<TrackingParticleCollection> trackingParticleHandle;
0247     event.getByLabel(edmCfg.getParameter<edm::InputTag>("trackingParticleTag"), trackingParticleHandle);
0248     LogTrace("l1tOmtfEventPrint") << "\nCandidateSimMuonMatcher::observeEventEnd trackingParticleHandle size "
0249                                   << trackingParticleHandle.product()->size() << std::endl;
0250 
0251     //TODO use other trackParticleFilter if needed <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
0252     std::function<bool(const TrackingParticle&)> trackParticleFilter =
0253         trackingParticleIsMuonInBx0;  //trackingParticleIsMuonInOmtfBx0;
0254     matchingResults =
0255         match(ghostBustedRegionalCands, ghostBustedProcMuons, trackingParticleHandle.product(), trackParticleFilter);
0256   }
0257 }
0258 
0259 void CandidateSimMuonMatcher::endJob() {}
0260 
0261 std::vector<const l1t::RegionalMuonCand*> CandidateSimMuonMatcher::ghostBust(
0262     const l1t::RegionalMuonCandBxCollection* mtfCands, const AlgoMuons& gbCandidates, AlgoMuons& ghostBustedProcMuons) {
0263   if (gbCandidates.size() != mtfCands->size(0)) {
0264     edm::LogError("l1tOmtfEventPrint") << "CandidateSimMuonMatcher::ghostBust(): gbCandidates.size() "
0265                                        << gbCandidates.size() << " != mtfCands.size() " << mtfCands->size();
0266   }
0267 
0268   boost::dynamic_bitset<> isKilled(mtfCands->size(0), false);
0269 
0270   for (unsigned int i1 = 0; i1 < mtfCands->size(0); ++i1) {
0271     if (mtfCands->at(0, i1).hwPt() == 0)
0272       continue;
0273     LogTrace("l1tOmtfEventPrint") << "\nCandidateSimMuonMatcher::ghostBust regionalCand pt " << std::setw(3)
0274                                   << mtfCands->at(0, i1).hwPt() << " qual " << std::setw(2)
0275                                   << mtfCands->at(0, i1).hwQual() << " proc " << std::setw(2)
0276                                   << mtfCands->at(0, i1).processor();
0277     for (unsigned int i2 = i1 + 1; i2 < mtfCands->size(0); ++i2) {
0278       auto& mtfCand1 = mtfCands->at(0, i1);
0279       auto& mtfCand2 = mtfCands->at(0, i2);
0280       if (mtfCand2.hwPt() == 0)
0281         continue;
0282 
0283       if (std::abs(mtfCand1.hwEta() - mtfCand2.hwEta()) < (0.3 / 0.010875)) {
0284         int gloablHwPhi1 = omtfConfig->calcGlobalPhi(mtfCand1.hwPhi(), mtfCand1.processor());
0285         int gloablHwPhi2 = omtfConfig->calcGlobalPhi(mtfCand2.hwPhi(), mtfCand2.processor());
0286 
0287         //folding phi
0288         int deltaPhi = std::abs(gloablHwPhi1 - gloablHwPhi2);
0289         if (deltaPhi > 576 / 2)
0290           deltaPhi = std::abs(deltaPhi - 576);
0291 
0292         //one can use the phi in radians like that:
0293         //double globalPhi1 = hwGmtPhiToGlobalPhi(omtfConfig->calcGlobalPhi( mtfCand1.hwPhi(), mtfCand1.processor() ) );
0294         //double globalPhi2 = hwGmtPhiToGlobalPhi(omtfConfig->calcGlobalPhi( mtfCand2.hwPhi(), mtfCand2.processor() ) );
0295 
0296         //0.0872664626 = 5 deg, i.e. the same window as in the OMTF ghost buster
0297         if (deltaPhi < 8) {
0298           //if (mtfCand1.hwQual() > mtfCand2.hwQual()) //TODO this is used in the uGMT
0299           //but this should be better - but probably the difference is not big
0300           if (gbCandidates[i1]->getFiredLayerCnt() > gbCandidates[i2]->getFiredLayerCnt()) {
0301             isKilled[i2] = true;
0302           } else
0303             isKilled[i1] = true;
0304         }
0305       }
0306     }
0307   }
0308 
0309   std::vector<const l1t::RegionalMuonCand*> resultCands;
0310 
0311   for (unsigned int i1 = 0; i1 < mtfCands->size(0); ++i1) {
0312     //dropping candidates with quality 0 !!!!!!!!!!!!!!!!!!!! fixme if not needed
0313     if (!isKilled[i1] && mtfCands->at(0, i1).hwPt()) {
0314       resultCands.push_back(&(mtfCands->at(0, i1)));
0315       ghostBustedProcMuons.push_back(gbCandidates.at(i1));
0316     }
0317 
0318     if (mtfCands->at(0, i1).hwPt()) {
0319       LogTrace("l1tOmtfEventPrint") << "CandidateSimMuonMatcher::ghostBust\n  regionalCand pt " << std::setw(3)
0320                                     << mtfCands->at(0, i1).hwPt() << " qual " << std::setw(2)
0321                                     << mtfCands->at(0, i1).hwQual() << " proc " << std::setw(2)
0322                                     << mtfCands->at(0, i1).processor() << " eta " << std::setw(4)
0323                                     << mtfCands->at(0, i1).hwEta() << " gloablEta " << std::setw(8)
0324                                     << mtfCands->at(0, i1).hwEta() * 0.010875 << " hwPhi " << std::setw(3)
0325                                     << mtfCands->at(0, i1).hwPhi() << " globalPhi " << std::setw(8)
0326                                     << hwGmtPhiToGlobalPhi(omtfConfig->calcGlobalPhi(mtfCands->at(0, i1).hwPhi(),
0327                                                                                      mtfCands->at(0, i1).processor()))
0328                                     << " fireadLayers " << std::bitset<18>(mtfCands->at(0, i1).trackAddress().at(0))
0329                                     << " gb isKilled " << isKilled.test(i1) << std::endl;
0330 
0331       LogTrace("l1tOmtfEventPrint") << *(gbCandidates.at(i1)) << std::endl;
0332     }
0333   }
0334 
0335   if (resultCands.size() >= 3)
0336     LogTrace("l1tOmtfEventPrint") << " ghost !!!!!! " << std::endl;
0337   LogTrace("l1tOmtfEventPrint") << std::endl;
0338 
0339   return resultCands;
0340 }
0341 
0342 TrajectoryStateOnSurface CandidateSimMuonMatcher::atStation2(const FreeTrajectoryState& ftsStart) const {
0343   // propagate to MB1, which defines the OMTF region (W+-2 MB1 is connected only to the OMTF)
0344   // 415 cm is R of RB1in, 660.5cm is |z| of the edge of MB2 (B field on)
0345   ReferenceCountingPointer<Surface> rpc = ReferenceCountingPointer<Surface>(
0346       new BoundCylinder(GlobalPoint(0., 0., 0.), TkRotation<float>(), SimpleCylinderBounds(415., 415., -660.5, 660.5)));
0347   TrajectoryStateOnSurface trackAtRPC = propagator->propagate(ftsStart, *rpc);
0348 
0349   return trackAtRPC;
0350 }
0351 
0352 FreeTrajectoryState CandidateSimMuonMatcher::simTrackToFts(const SimTrack& simTrackPtr, const SimVertex& simVertex) {
0353   int charge = simTrackPtr.type() > 0 ? -1 : 1;  //works for muons
0354 
0355   GlobalVector p3GV(simTrackPtr.momentum().x(), simTrackPtr.momentum().y(), simTrackPtr.momentum().z());
0356   GlobalPoint r3GP(simVertex.position().x(), simVertex.position().y(), simVertex.position().z());
0357 
0358   GlobalTrajectoryParameters tPars(r3GP, p3GV, charge, &*magField);
0359 
0360   return FreeTrajectoryState(tPars);
0361 }
0362 
0363 FreeTrajectoryState CandidateSimMuonMatcher::simTrackToFts(const TrackingParticle& trackingParticle) {
0364   int charge = trackingParticle.pdgId() > 0 ? -1 : 1;  //works for muons
0365 
0366   GlobalVector p3GV(trackingParticle.momentum().x(), trackingParticle.momentum().y(), trackingParticle.momentum().z());
0367   GlobalPoint r3GP(trackingParticle.vx(), trackingParticle.vy(), trackingParticle.vz());
0368 
0369   GlobalTrajectoryParameters tPars(r3GP, p3GV, charge, &*magField);
0370 
0371   return FreeTrajectoryState(tPars);
0372 }
0373 
0374 TrajectoryStateOnSurface CandidateSimMuonMatcher::propagate(const SimTrack& simTrack,
0375                                                             const edm::SimVertexContainer* simVertices) {
0376   SimVertex simVertex;
0377   int vtxInd = simTrack.vertIndex();
0378   if (vtxInd < 0) {
0379     edm::LogImportant("l1tOmtfEventPrint") << "Track with no vertex, defaulting to (0,0,0)";
0380   } else {
0381     simVertex = simVertices->at(vtxInd);
0382     if (((int)simVertex.vertexId()) != vtxInd) {
0383       edm::LogImportant("l1tOmtfEventPrint") << "simVertex.vertexId() != vtxInd. simVertex.vertexId() "
0384                                              << simVertex.vertexId() << " vtxInd " << vtxInd << " !!!!!!!!!!!!!!!!!";
0385     }
0386   }
0387 
0388   FreeTrajectoryState ftsTrack = simTrackToFts(simTrack, simVertex);
0389 
0390   TrajectoryStateOnSurface tsof = atStation2(ftsTrack);  //propagation
0391 
0392   return tsof;
0393 }
0394 
0395 TrajectoryStateOnSurface CandidateSimMuonMatcher::propagate(const TrackingParticle& trackingParticle) {
0396   FreeTrajectoryState ftsTrack = simTrackToFts(trackingParticle);
0397 
0398   TrajectoryStateOnSurface tsof = atStation2(ftsTrack);  //propagation
0399 
0400   return tsof;
0401 }
0402 
0403 float normal_pdf(float x, float m, float s) {
0404   static const float inv_sqrt_2pi = 0.3989422804014327;
0405   float a = (x - m) / s;
0406 
0407   return inv_sqrt_2pi / s * std::exp(-0.5 * a * a);
0408 }
0409 
0410 MatchingResult CandidateSimMuonMatcher::match(const l1t::RegionalMuonCand* muonCand,
0411                                               const AlgoMuonPtr& procMuon,
0412                                               const SimTrack& simTrack,
0413                                               TrajectoryStateOnSurface& tsof) {
0414   MatchingResult result(simTrack);
0415 
0416   double candGloablEta = muonCand->hwEta() * 0.010875;
0417   //if (std::abs(simTrack.momentum().eta() - candGloablEta) < 0.3) //has no sense for displaced muons
0418   {
0419     double candGlobalPhi = omtfConfig->calcGlobalPhi(muonCand->hwPhi(), muonCand->processor());
0420     candGlobalPhi = hwGmtPhiToGlobalPhi(candGlobalPhi);
0421 
0422     if (candGlobalPhi > M_PI)
0423       candGlobalPhi = candGlobalPhi - (2. * M_PI);
0424 
0425     result.deltaPhi = foldPhi(tsof.globalPosition().phi() - candGlobalPhi);
0426     result.deltaEta = tsof.globalPosition().eta() - candGloablEta;
0427 
0428     result.propagatedPhi = tsof.globalPosition().phi();
0429     result.propagatedEta = tsof.globalPosition().eta();
0430 
0431     double mean = 0;
0432     double sigma = 1;
0433     //if(!fillMean)
0434     {
0435       auto ptBin = deltaPhiPropCandMean->FindBin(simTrack.momentum().pt());
0436       mean = deltaPhiPropCandMean->GetBinContent(ptBin);
0437       sigma = deltaPhiPropCandStdDev->GetBinContent(ptBin);
0438     }
0439     result.matchingLikelihood = normal_pdf(result.deltaPhi, mean, sigma);  //TODO temporary solution
0440 
0441     result.muonCand = muonCand;
0442     result.procMuon = procMuon;
0443 
0444     double treshold = 6. * sigma;
0445     if (simTrack.momentum().pt() > 20)
0446       treshold = 7. * sigma;
0447     if (simTrack.momentum().pt() > 100)
0448       treshold = 20. * sigma;
0449 
0450     //for displaced muons in H2ll
0451     treshold = 0.15;  //pt > 30
0452     if (simTrack.momentum().pt() <
0453         10)  //TODO!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! tune the threshold!!!!!!
0454       treshold = 0.3;
0455     else if (simTrack.momentum().pt() <
0456              30)  //TODO!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! tune the threshold!!!!!!
0457       treshold = 0.22;
0458 
0459     if (std::abs(result.deltaPhi - mean) < treshold && std::abs(result.deltaEta) < 0.3)
0460       result.result = MatchingResult::ResultType::matched;
0461 
0462     LogTrace("l1tOmtfEventPrint") << "CandidateSimMuonMatcher::match: simTrack type " << simTrack.type() << " pt "
0463                                   << std::setw(8) << simTrack.momentum().pt() << " eta " << std::setw(8)
0464                                   << simTrack.momentum().eta() << " phi " << std::setw(8) << simTrack.momentum().phi()
0465                                   << " propagation eta " << std::setw(8) << tsof.globalPosition().eta() << " phi "
0466                                   << tsof.globalPosition().phi() << "\n             muonCand pt " << std::setw(8)
0467                                   << muonCand->hwPt() << " candGloablEta " << std::setw(8) << candGloablEta
0468                                   << " candGlobalPhi " << std::setw(8) << candGlobalPhi << " hwQual "
0469                                   << muonCand->hwQual() << " deltaEta " << std::setw(8) << result.deltaEta
0470                                   << " deltaPhi " << std::setw(8) << result.deltaPhi << " sigma " << std::setw(8)
0471                                   << sigma << " Likelihood " << std::setw(8) << result.matchingLikelihood << " result "
0472                                   << (short)result.result << std::endl;
0473   }
0474 
0475   return result;
0476 }
0477 
0478 MatchingResult CandidateSimMuonMatcher::match(const l1t::RegionalMuonCand* muonCand,
0479                                               const AlgoMuonPtr& procMuon,
0480                                               const TrackingParticle& trackingParticle,
0481                                               TrajectoryStateOnSurface& tsof) {
0482   MatchingResult result(trackingParticle);
0483 
0484   double candGloablEta = muonCand->hwEta() * 0.010875;
0485   //if (std::abs(trackingParticle.momentum().eta() - candGloablEta) < 0.3)  //has no sense for displaced muons
0486   {
0487     double candGlobalPhi = omtfConfig->calcGlobalPhi(muonCand->hwPhi(), muonCand->processor());
0488     candGlobalPhi = hwGmtPhiToGlobalPhi(candGlobalPhi);
0489 
0490     if (candGlobalPhi > M_PI)
0491       candGlobalPhi = candGlobalPhi - (2. * M_PI);
0492 
0493     result.deltaPhi = foldPhi(tsof.globalPosition().phi() - candGlobalPhi);
0494     result.deltaEta = tsof.globalPosition().eta() - candGloablEta;
0495 
0496     result.propagatedPhi = tsof.globalPosition().phi();
0497     result.propagatedEta = tsof.globalPosition().eta();
0498 
0499     double mean = 0;
0500     double sigma = 1;
0501     //if(!fillMean)
0502     {
0503       auto ptBin = deltaPhiPropCandMean->FindBin(trackingParticle.pt());
0504 
0505       mean = deltaPhiPropCandMean->GetBinContent(ptBin);
0506       sigma = deltaPhiPropCandStdDev->GetBinContent(ptBin);
0507     }
0508 
0509     result.matchingLikelihood = normal_pdf(result.deltaPhi, mean, sigma);  //TODO temporary solution
0510 
0511     result.muonCand = muonCand;
0512     result.procMuon = procMuon;
0513 
0514     double treshold = 6. * sigma;
0515     if (trackingParticle.pt() > 20)
0516       treshold = 7. * sigma;
0517     if (trackingParticle.pt() > 100)
0518       treshold = 20. * sigma;
0519 
0520     if (std::abs(result.deltaPhi - mean) < treshold && std::abs(result.deltaEta) < 0.3)
0521       result.result = MatchingResult::ResultType::matched;
0522 
0523     LogTrace("l1tOmtfEventPrint") << "CandidateSimMuonMatcher::match: trackingParticle type "
0524                                   << trackingParticle.pdgId() << " pt " << std::setw(8) << trackingParticle.pt()
0525                                   << " eta " << std::setw(8) << trackingParticle.momentum().eta() << " phi "
0526                                   << std::setw(8) << trackingParticle.momentum().phi() << " propagation eta "
0527                                   << std::setw(8) << tsof.globalPosition().eta() << " phi "
0528                                   << tsof.globalPosition().phi() << " muonCand pt " << std::setw(8) << muonCand->hwPt()
0529                                   << " candGloablEta " << std::setw(8) << candGloablEta << " candGlobalPhi "
0530                                   << std::setw(8) << candGlobalPhi << " hwQual " << muonCand->hwQual() << " deltaEta "
0531                                   << std::setw(8) << result.deltaEta << " deltaPhi " << std::setw(8) << result.deltaPhi
0532                                   << " Likelihood " << std::setw(8) << result.matchingLikelihood << " result "
0533                                   << (short)result.result << std::endl;
0534   }
0535 
0536   return result;
0537 }
0538 
0539 std::vector<MatchingResult> CandidateSimMuonMatcher::cleanMatching(std::vector<MatchingResult> matchingResults,
0540                                                                    std::vector<const l1t::RegionalMuonCand*>& muonCands,
0541                                                                    AlgoMuons& ghostBustedProcMuons) {
0542   //Cleaning the matching
0543   std::sort(
0544       matchingResults.begin(), matchingResults.end(), [](const MatchingResult& a, const MatchingResult& b) -> bool {
0545         return a.matchingLikelihood > b.matchingLikelihood;
0546       });
0547 
0548   for (unsigned int i1 = 0; i1 < matchingResults.size(); i1++) {
0549     if (matchingResults[i1].result == MatchingResult::ResultType::matched) {
0550       for (unsigned int i2 = i1 + 1; i2 < matchingResults.size(); i2++) {
0551         if ((matchingResults[i1].trackingParticle &&
0552              matchingResults[i1].trackingParticle == matchingResults[i2].trackingParticle) ||
0553             (matchingResults[i1].simTrack && matchingResults[i1].simTrack == matchingResults[i2].simTrack) ||
0554             (matchingResults[i1].muonCand == matchingResults[i2].muonCand)) {
0555           //if matchingResults[i1].muonCand == false, then it is also OK here
0556           matchingResults[i2].result = MatchingResult::ResultType::duplicate;
0557         }
0558       }
0559     }
0560   }
0561 
0562   std::vector<MatchingResult> cleanedMatchingResults;
0563   for (auto& matchingResult : matchingResults) {
0564     if (matchingResult.result == MatchingResult::ResultType::matched || matchingResult.muonCand == nullptr)
0565       //adding also the simTracks that are not matched at all, before it is assured that they are not duplicates
0566       cleanedMatchingResults.push_back(matchingResult);
0567     if (matchingResult.result == MatchingResult::ResultType::matched) {
0568       /* TODO fix if filling of the deltaPhiPropCandMean and deltaPhiPropCandStdDev is needed
0569       if(fillMean) {
0570       double ptGen  = matchingResult.genPt;
0571         deltaPhiPropCandMean->Fill(ptGen, matchingResult.deltaPhi); //filling overflow is ok here
0572         deltaPhiPropCandStdDev->Fill(ptGen, matchingResult.deltaPhi * matchingResult.deltaPhi);
0573       }*/
0574     }
0575   }
0576 
0577   //adding the muonCand-s that were not matched, i.e. in order to analyze them later
0578   unsigned int iCand = 0;
0579   for (auto& muonCand : muonCands) {
0580     bool isMatched = false;
0581     for (auto& matchingResult : cleanedMatchingResults) {
0582       if (matchingResult.muonCand == muonCand) {
0583         isMatched = true;
0584         break;
0585       }
0586     }
0587 
0588     if (!isMatched) {
0589       MatchingResult result;
0590       result.muonCand = muonCand;
0591       result.procMuon = ghostBustedProcMuons.at(iCand);
0592       cleanedMatchingResults.push_back(result);
0593     }
0594     iCand++;
0595   }
0596 
0597   LogTrace("l1tOmtfEventPrint") << "CandidateSimMuonMatcher::cleanMatching:" << __LINE__
0598                                 << " cleanedMatchingResults:" << std::endl;
0599   for (auto& result : cleanedMatchingResults) {
0600     if (result.trackingParticle || result.simTrack)
0601       LogTrace("l1tOmtfEventPrint") << "CandidateSimMuonMatcher::cleanMatching:" << __LINE__ << " simTrack type "
0602                                     << result.pdgId << " pt " << std::setw(8) << result.genPt << " eta " << std::setw(8)
0603                                     << result.genEta << " phi " << std::setw(8) << result.genPhi;
0604     else
0605       LogTrace("l1tOmtfEventPrint") << "no matched track ";
0606 
0607     if (result.muonCand) {
0608       LogTrace("l1tOmtfEventPrint") << " muonCand pt " << std::setw(8) << result.muonCand->hwPt() << " hwQual "
0609                                     << result.muonCand->hwQual() << " hwEta " << result.muonCand->hwEta()
0610                                     << " deltaEta " << std::setw(8) << result.deltaEta << " deltaPhi " << std::setw(8)
0611                                     << result.deltaPhi << " Likelihood " << std::setw(8) << result.matchingLikelihood
0612                                     << " result " << (short)result.result;
0613       LogTrace("l1tOmtfEventPrint") << " procMuon " << *(result.procMuon) << std::endl;
0614     } else
0615       LogTrace("l1tOmtfEventPrint") << " no muonCand "
0616                                     << " result " << (short)result.result << std::endl;
0617   }
0618   LogTrace("l1tOmtfEventPrint") << " " << std::endl;
0619 
0620   return cleanedMatchingResults;
0621 }
0622 
0623 std::vector<MatchingResult> CandidateSimMuonMatcher::match(std::vector<const l1t::RegionalMuonCand*>& muonCands,
0624                                                            AlgoMuons& ghostBustedProcMuons,
0625                                                            const edm::SimTrackContainer* simTracks,
0626                                                            const edm::SimVertexContainer* simVertices,
0627                                                            std::function<bool(const SimTrack&)> const& simTrackFilter) {
0628   std::vector<MatchingResult> matchingResults;
0629 
0630   for (auto& simTrack : *simTracks) {
0631     if (!simTrackFilter(simTrack))
0632       continue;
0633 
0634     LogTrace("l1tOmtfEventPrint") << "CandidateSimMuonMatcher::match, simTrack type " << std::setw(3) << simTrack.type()
0635                                   << " pt " << std::setw(9) << simTrack.momentum().pt() << " eta " << std::setw(9)
0636                                   << simTrack.momentum().eta() << " phi " << std::setw(9) << simTrack.momentum().phi()
0637                                   << std::endl;
0638 
0639     bool matched = false;
0640 
0641     TrajectoryStateOnSurface tsof = propagate(simTrack, simVertices);
0642     if (!tsof.isValid()) {  //no sense to do matching
0643       MatchingResult result(simTrack);
0644       result.result = MatchingResult::ResultType::propagationFailed;
0645       LogTrace("l1tOmtfEventPrint") << __FUNCTION__ << ":" << __LINE__ << " propagation failed: genPt " << result.genPt
0646                                     << " genEta " << result.genEta << " eventId " << simTrack.eventId().event()
0647                                     << std::endl;
0648 
0649       matchingResults.push_back(result);
0650     } else {
0651       //checking if the propagated track is inside the OMTF range, TODO - tune the range!!!!!!!!!!!!!!!!!
0652       //eta 0.7 is the beginning of the MB2,
0653       //the eta range wider than the nominal OMTF region is needed, as in any case muons outside this region are seen by the OMTF
0654       //so it better to train the nn suich that is able to measure its pt, as it may affect the rate
0655       if ((std::abs(tsof.globalPosition().eta()) >= 0.7) && (std::abs(tsof.globalPosition().eta()) <= 1.31)) {
0656         LogTrace("l1tOmtfEventPrint")
0657             << "CandidateSimMuonMatcher::match simTrack IS in OMTF region, matching to the omtfCands";
0658       } else {
0659         LogTrace("l1tOmtfEventPrint") << "simTrack NOT in OMTF region ";
0660         continue;
0661       }
0662 
0663       /* TODO fix if filling of the deltaPhiPropCandMean and deltaPhiPropCandStdDev is needed
0664     double ptGen = simTrack.momentum().pt();
0665     if(ptGen >= deltaPhiVertexProp->GetXaxis()->GetXmax())
0666       ptGen = deltaPhiVertexProp->GetXaxis()->GetXmax() - 0.01;
0667 
0668     deltaPhiVertexProp->Fill(ptGen, simTrack.momentum().phi() - tsof.globalPosition().phi());*/
0669 
0670       unsigned int iCand = 0;
0671       for (auto& muonCand : muonCands) {
0672         //dropping very low quality candidates, as they are fakes usually - but it has no sense, then the results are not conclusive
0673         //if(muonCand->hwQual() > 1)
0674         {
0675           MatchingResult result;
0676           if (tsof.isValid()) {
0677             result = match(muonCand, ghostBustedProcMuons.at(iCand), simTrack, tsof);
0678           }
0679           int vtxInd = simTrack.vertIndex();
0680           if (vtxInd >= 0) {
0681             result.simVertex = &(simVertices->at(
0682                 vtxInd));  //TODO ?????? something strange is here, was commented in the previous version
0683           }
0684           if (result.result == MatchingResult::ResultType::matched) {
0685             matchingResults.push_back(result);
0686             matched = true;
0687           }
0688         }
0689         iCand++;
0690       }
0691 
0692       if (!matched) {  //we are adding also if it was not matching to any candidate
0693         MatchingResult result(simTrack);
0694         matchingResults.push_back(result);
0695         LogTrace("l1tOmtfEventPrint") << __FUNCTION__ << ":" << __LINE__ << " no matching candidate found" << std::endl;
0696       }
0697     }
0698   }
0699 
0700   return cleanMatching(matchingResults, muonCands, ghostBustedProcMuons);
0701 }
0702 
0703 std::vector<MatchingResult> CandidateSimMuonMatcher::match(
0704     std::vector<const l1t::RegionalMuonCand*>& muonCands,
0705     AlgoMuons& ghostBustedProcMuons,
0706     const TrackingParticleCollection* trackingParticles,
0707     std::function<bool(const TrackingParticle&)> const& simTrackFilter) {
0708   std::vector<MatchingResult> matchingResults;
0709   LogTrace("l1tOmtfEventPrint") << "CandidateSimMuonMatcher::match trackingParticles->size() "
0710                                 << trackingParticles->size() << std::endl;
0711 
0712   for (auto& trackingParticle : *trackingParticles) {
0713     //LogTrace("l1tOmtfEventPrint") <<"CandidateSimMuonMatcher::match:"<<__LINE__<<" trackingParticle type "<<std::setw(3)<<trackingParticle.pdgId()<<" pt "<<std::setw(9)<<trackingParticle.pt()<<" eta "<<std::setw(9)<<trackingParticle.momentum().eta()<<" phi "<<std::setw(9)<<trackingParticle.momentum().phi()<<std::endl;
0714 
0715     if (simTrackFilter(trackingParticle) == false)
0716       continue;
0717 
0718     LogTrace("l1tOmtfEventPrint") << "CandidateSimMuonMatcher::match, trackingParticle type " << std::setw(3)
0719                                   << trackingParticle.pdgId() << " pt " << std::setw(9) << trackingParticle.pt()
0720                                   << " eta " << std::setw(9) << trackingParticle.momentum().eta() << " phi "
0721                                   << std::setw(9) << trackingParticle.momentum().phi() << std::endl;
0722 
0723     bool matched = false;
0724 
0725     TrajectoryStateOnSurface tsof = propagate(trackingParticle);
0726     if (!tsof.isValid()) {
0727       LogTrace("l1tOmtfEventPrint") << "CandidateSimMuonMatcher::match:" << __LINE__ << " propagation failed"
0728                                     << std::endl;
0729       MatchingResult result;
0730       result.result = MatchingResult::ResultType::propagationFailed;
0731       continue;  //no sense to do matching
0732     }
0733 
0734     LogTrace("l1tOmtfEventPrint") << "CandidateSimMuonMatcher::match, tsof.globalPosition().eta() "
0735                                   << tsof.globalPosition().eta();
0736 
0737     //checking if the propagated track is inside the OMTF range, TODO - tune the range!!!!!!!!!!!!!!!!!
0738     //eta 0.7 is the beginning of the MB2, while 1.31 is mid of RE2 + some margin
0739     if ((std::abs(tsof.globalPosition().eta()) >= 0.7) && (std::abs(tsof.globalPosition().eta()) <= 1.31)) {
0740       LogTrace("l1tOmtfEventPrint")
0741           << "CandidateSimMuonMatcher::match trackingParticle IS in OMTF region, matching to the omtfCands";
0742     } else {
0743       LogTrace("l1tOmtfEventPrint") << "trackingParticle NOT in OMTF region ";
0744       continue;
0745     }
0746 
0747     /* TODO fix if filling of the deltaPhiPropCandMean and deltaPhiPropCandStdDev is needed
0748     double ptGen = trackingParticle.pt();
0749     if(ptGen >= deltaPhiVertexProp->GetXaxis()->GetXmax())
0750       ptGen = deltaPhiVertexProp->GetXaxis()->GetXmax() - 0.01;
0751 
0752     deltaPhiVertexProp->Fill(ptGen, trackingParticle.momentum().phi() - tsof.globalPosition().phi());
0753      */
0754 
0755     unsigned int iCand = 0;
0756     for (auto& muonCand : muonCands) {
0757       //dropping very low quality candidates, as they are fakes usually - but it has no sense, then the results are not conclusive then
0758       /*if(muonCand->hwQual() <= 1)
0759         continue; */
0760 
0761       MatchingResult result;
0762       if (tsof.isValid()) {
0763         result = match(muonCand, ghostBustedProcMuons.at(iCand), trackingParticle, tsof);
0764       }
0765       iCand++;
0766 
0767       if (result.result == MatchingResult::ResultType::matched) {
0768         matchingResults.push_back(result);
0769         matched = true;
0770       }
0771     }
0772 
0773     if (!matched) {  //we are adding result also if it there was no matching to any candidate
0774       MatchingResult result(trackingParticle);
0775       matchingResults.push_back(result);
0776       LogTrace("l1tOmtfEventPrint") << __FUNCTION__ << ":" << __LINE__ << " no matching candidate found" << std::endl;
0777     }
0778   }
0779 
0780   return cleanMatching(matchingResults, muonCands, ghostBustedProcMuons);
0781 }
0782 
0783 std::vector<MatchingResult> CandidateSimMuonMatcher::matchSimple(
0784     std::vector<const l1t::RegionalMuonCand*>& muonCands,
0785     AlgoMuons& ghostBustedProcMuons,
0786     const edm::SimTrackContainer* simTracks,
0787     const edm::SimVertexContainer* simVertices,
0788     std::function<bool(const SimTrack&)> const& simTrackFilter) {
0789   std::vector<MatchingResult> matchingResults;
0790 
0791   for (auto& simTrack : *simTracks) {
0792     if (!simTrackFilter(simTrack))
0793       continue;
0794 
0795     LogTrace("l1tOmtfEventPrint") << "CandidateSimMuonMatcher::matchSimple: simTrack type " << std::setw(3)
0796                                   << simTrack.type() << " pt " << std::setw(9) << simTrack.momentum().pt() << " eta "
0797                                   << std::setw(9) << simTrack.momentum().eta() << " phi " << std::setw(9)
0798                                   << simTrack.momentum().phi() << std::endl;
0799 
0800     bool matched = false;
0801 
0802     unsigned int iCand = 0;
0803     for (auto& muonCand : muonCands) {
0804       //dropping very low quality candidates, as they are fakes usually - but it has no sense, then the results are not conclusive
0805       //if(muonCand->hwQual() > 1)
0806       {
0807         MatchingResult result(simTrack);
0808 
0809         double candGloablEta = muonCand->hwEta() * 0.010875;
0810         double candGlobalPhi = omtfConfig->calcGlobalPhi(muonCand->hwPhi(), muonCand->processor());
0811         candGlobalPhi = hwGmtPhiToGlobalPhi(candGlobalPhi);
0812 
0813         if (candGlobalPhi > M_PI)
0814           candGlobalPhi = candGlobalPhi - (2. * M_PI);
0815 
0816         result.deltaPhi = foldPhi(result.genPhi - candGlobalPhi);
0817         result.deltaEta = result.genEta - candGloablEta;
0818 
0819         result.propagatedPhi = 0;
0820         result.propagatedEta = 0;
0821 
0822         result.muonCand = muonCand;
0823         result.procMuon = ghostBustedProcMuons.at(iCand);
0824 
0825         //TODO histogram can be used, like in the usercode/L1MuonAnalyzer  MuonMatcher::matchWithoutPorpagation
0826         //for prompt muons
0827         /*double treshold = 0.3;
0828         if (simTrack.momentum().pt() < 5)
0829           treshold = 1.5;
0830         else if (simTrack.momentum().pt() < 8)
0831           treshold = 0.8;
0832         else if (simTrack.momentum().pt() < 10)
0833           treshold = 0.6;
0834         else if (simTrack.momentum().pt() < 20)
0835           treshold = 0.5;*/
0836 
0837         //for displaced muons
0838         double treshold = 0.7;
0839         if (simTrack.momentum().pt() < 5)
0840           treshold = 1.5;
0841         else if (simTrack.momentum().pt() < 10)
0842           treshold = 1.0;
0843         else if (simTrack.momentum().pt() < 20)
0844           treshold = 0.7;
0845 
0846         if (std::abs(result.deltaPhi) < treshold && std::abs(result.deltaEta) < 0.5) {
0847           result.result = MatchingResult::ResultType::matched;
0848           //matchingLikelihood is needed in the cleanMatching, so we put something
0849           if (std::abs(result.deltaPhi) < 0.001)
0850             result.matchingLikelihood = 1. / 0.001;
0851           else
0852             result.matchingLikelihood = 1. / std::abs(result.deltaPhi);
0853         }
0854 
0855         LogTrace("l1tOmtfEventPrint") << "CandidateSimMuonMatcher::matchSimple: simTrack type " << simTrack.type()
0856                                       << " pt " << std::setw(8) << simTrack.momentum().pt() << " eta " << std::setw(8)
0857                                       << simTrack.momentum().eta() << " phi " << std::setw(8)
0858                                       << simTrack.momentum().phi() << "\n             muonCand pt " << std::setw(8)
0859                                       << muonCand->hwPt() << " candGloablEta " << std::setw(8) << candGloablEta
0860                                       << " candGlobalPhi " << std::setw(8) << candGlobalPhi << " hwQual "
0861                                       << muonCand->hwQual() << " deltaEta " << std::setw(8) << result.deltaEta
0862                                       << " deltaPhi " << std::setw(8) << result.deltaPhi << " matchingLikelihood "
0863                                       << result.matchingLikelihood << " result " << (short)result.result << std::endl;
0864 
0865         int vtxInd = simTrack.vertIndex();
0866         if (vtxInd >= 0) {
0867           result.simVertex = &(
0868               simVertices->at(vtxInd));  //TODO ?????? something strange is here, was commented in the previous version
0869         }
0870         if (result.result == MatchingResult::ResultType::matched) {
0871           matchingResults.push_back(result);
0872           matched = true;
0873         }
0874       }
0875       iCand++;
0876     }
0877 
0878     if (!matched) {  //we are adding also if it was not matched to any candidate
0879       MatchingResult result(simTrack);
0880       matchingResults.push_back(result);
0881       LogTrace("l1tOmtfEventPrint") << __FUNCTION__ << ":" << __LINE__ << " no matching candidate found" << std::endl;
0882     }
0883   }
0884 
0885   return cleanMatching(matchingResults, muonCands, ghostBustedProcMuons);
0886 }