Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-10-08 05:11:49

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     //for displaced muons in H2ll
0445     double treshold = 0.15;  //pt > 30
0446     if (simTrack.momentum().pt() <
0447         10)  //TODO!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! tune the threshold!!!!!!
0448       treshold = 0.3;
0449     else if (simTrack.momentum().pt() <
0450              30)  //TODO!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! tune the threshold!!!!!!
0451       treshold = 0.22;
0452 
0453     if (std::abs(result.deltaPhi - mean) < treshold && std::abs(result.deltaEta) < 0.3)
0454       result.result = MatchingResult::ResultType::matched;
0455 
0456     LogTrace("l1tOmtfEventPrint") << "CandidateSimMuonMatcher::match: simTrack type " << simTrack.type() << " pt "
0457                                   << std::setw(8) << simTrack.momentum().pt() << " eta " << std::setw(8)
0458                                   << simTrack.momentum().eta() << " phi " << std::setw(8) << simTrack.momentum().phi()
0459                                   << " propagation eta " << std::setw(8) << tsof.globalPosition().eta() << " phi "
0460                                   << tsof.globalPosition().phi() << "\n             muonCand pt " << std::setw(8)
0461                                   << muonCand->hwPt() << " candGloablEta " << std::setw(8) << candGloablEta
0462                                   << " candGlobalPhi " << std::setw(8) << candGlobalPhi << " hwQual "
0463                                   << muonCand->hwQual() << " deltaEta " << std::setw(8) << result.deltaEta
0464                                   << " deltaPhi " << std::setw(8) << result.deltaPhi << " sigma " << std::setw(8)
0465                                   << sigma << " Likelihood " << std::setw(8) << result.matchingLikelihood << " result "
0466                                   << (short)result.result << std::endl;
0467   }
0468 
0469   return result;
0470 }
0471 
0472 MatchingResult CandidateSimMuonMatcher::match(const l1t::RegionalMuonCand* muonCand,
0473                                               const AlgoMuonPtr& procMuon,
0474                                               const TrackingParticle& trackingParticle,
0475                                               TrajectoryStateOnSurface& tsof) {
0476   MatchingResult result(trackingParticle);
0477 
0478   double candGloablEta = muonCand->hwEta() * 0.010875;
0479   //if (std::abs(trackingParticle.momentum().eta() - candGloablEta) < 0.3)  //has no sense for displaced muons
0480   {
0481     double candGlobalPhi = omtfConfig->calcGlobalPhi(muonCand->hwPhi(), muonCand->processor());
0482     candGlobalPhi = hwGmtPhiToGlobalPhi(candGlobalPhi);
0483 
0484     if (candGlobalPhi > M_PI)
0485       candGlobalPhi = candGlobalPhi - (2. * M_PI);
0486 
0487     result.deltaPhi = foldPhi(tsof.globalPosition().phi() - candGlobalPhi);
0488     result.deltaEta = tsof.globalPosition().eta() - candGloablEta;
0489 
0490     result.propagatedPhi = tsof.globalPosition().phi();
0491     result.propagatedEta = tsof.globalPosition().eta();
0492 
0493     double mean = 0;
0494     double sigma = 1;
0495     //if(!fillMean)
0496     {
0497       auto ptBin = deltaPhiPropCandMean->FindBin(trackingParticle.pt());
0498 
0499       mean = deltaPhiPropCandMean->GetBinContent(ptBin);
0500       sigma = deltaPhiPropCandStdDev->GetBinContent(ptBin);
0501     }
0502 
0503     result.matchingLikelihood = normal_pdf(result.deltaPhi, mean, sigma);  //TODO temporary solution
0504 
0505     result.muonCand = muonCand;
0506     result.procMuon = procMuon;
0507 
0508     double treshold = 6. * sigma;
0509     if (trackingParticle.pt() > 20)
0510       treshold = 7. * sigma;
0511     if (trackingParticle.pt() > 100)
0512       treshold = 20. * sigma;
0513 
0514     if (std::abs(result.deltaPhi - mean) < treshold && std::abs(result.deltaEta) < 0.3)
0515       result.result = MatchingResult::ResultType::matched;
0516 
0517     LogTrace("l1tOmtfEventPrint") << "CandidateSimMuonMatcher::match: trackingParticle type "
0518                                   << trackingParticle.pdgId() << " pt " << std::setw(8) << trackingParticle.pt()
0519                                   << " eta " << std::setw(8) << trackingParticle.momentum().eta() << " phi "
0520                                   << std::setw(8) << trackingParticle.momentum().phi() << " propagation eta "
0521                                   << std::setw(8) << tsof.globalPosition().eta() << " phi "
0522                                   << tsof.globalPosition().phi() << " muonCand pt " << std::setw(8) << muonCand->hwPt()
0523                                   << " candGloablEta " << std::setw(8) << candGloablEta << " candGlobalPhi "
0524                                   << std::setw(8) << candGlobalPhi << " hwQual " << muonCand->hwQual() << " deltaEta "
0525                                   << std::setw(8) << result.deltaEta << " deltaPhi " << std::setw(8) << result.deltaPhi
0526                                   << " Likelihood " << std::setw(8) << result.matchingLikelihood << " result "
0527                                   << (short)result.result << std::endl;
0528   }
0529 
0530   return result;
0531 }
0532 
0533 std::vector<MatchingResult> CandidateSimMuonMatcher::cleanMatching(std::vector<MatchingResult> matchingResults,
0534                                                                    std::vector<const l1t::RegionalMuonCand*>& muonCands,
0535                                                                    AlgoMuons& ghostBustedProcMuons) {
0536   //Cleaning the matching
0537   std::sort(
0538       matchingResults.begin(), matchingResults.end(), [](const MatchingResult& a, const MatchingResult& b) -> bool {
0539         return a.matchingLikelihood > b.matchingLikelihood;
0540       });
0541 
0542   for (unsigned int i1 = 0; i1 < matchingResults.size(); i1++) {
0543     if (matchingResults[i1].result == MatchingResult::ResultType::matched) {
0544       for (unsigned int i2 = i1 + 1; i2 < matchingResults.size(); i2++) {
0545         if ((matchingResults[i1].trackingParticle &&
0546              matchingResults[i1].trackingParticle == matchingResults[i2].trackingParticle) ||
0547             (matchingResults[i1].simTrack && matchingResults[i1].simTrack == matchingResults[i2].simTrack) ||
0548             (matchingResults[i1].muonCand == matchingResults[i2].muonCand)) {
0549           //if matchingResults[i1].muonCand == false, then it is also OK here
0550           matchingResults[i2].result = MatchingResult::ResultType::duplicate;
0551         }
0552       }
0553     }
0554   }
0555 
0556   std::vector<MatchingResult> cleanedMatchingResults;
0557   for (auto& matchingResult : matchingResults) {
0558     if (matchingResult.result == MatchingResult::ResultType::matched || matchingResult.muonCand == nullptr)
0559       //adding also the simTracks that are not matched at all, before it is assured that they are not duplicates
0560       cleanedMatchingResults.push_back(matchingResult);
0561     if (matchingResult.result == MatchingResult::ResultType::matched) {
0562       /* TODO fix if filling of the deltaPhiPropCandMean and deltaPhiPropCandStdDev is needed
0563       if(fillMean) {
0564       double ptGen  = matchingResult.genPt;
0565         deltaPhiPropCandMean->Fill(ptGen, matchingResult.deltaPhi); //filling overflow is ok here
0566         deltaPhiPropCandStdDev->Fill(ptGen, matchingResult.deltaPhi * matchingResult.deltaPhi);
0567       }*/
0568     }
0569   }
0570 
0571   //adding the muonCand-s that were not matched, i.e. in order to analyze them later
0572   unsigned int iCand = 0;
0573   for (auto& muonCand : muonCands) {
0574     bool isMatched = false;
0575     for (auto& matchingResult : cleanedMatchingResults) {
0576       if (matchingResult.muonCand == muonCand) {
0577         isMatched = true;
0578         break;
0579       }
0580     }
0581 
0582     if (!isMatched) {
0583       MatchingResult result;
0584       result.muonCand = muonCand;
0585       result.procMuon = ghostBustedProcMuons.at(iCand);
0586       cleanedMatchingResults.push_back(result);
0587     }
0588     iCand++;
0589   }
0590 
0591   LogTrace("l1tOmtfEventPrint") << "CandidateSimMuonMatcher::cleanMatching:" << __LINE__
0592                                 << " cleanedMatchingResults:" << std::endl;
0593   for (auto& result : cleanedMatchingResults) {
0594     if (result.trackingParticle || result.simTrack)
0595       LogTrace("l1tOmtfEventPrint") << "CandidateSimMuonMatcher::cleanMatching:" << __LINE__ << " simTrack type "
0596                                     << result.pdgId << " pt " << std::setw(8) << result.genPt << " eta " << std::setw(8)
0597                                     << result.genEta << " phi " << std::setw(8) << result.genPhi;
0598     else
0599       LogTrace("l1tOmtfEventPrint") << "no matched track ";
0600 
0601     if (result.muonCand) {
0602       LogTrace("l1tOmtfEventPrint") << " muonCand pt " << std::setw(8) << result.muonCand->hwPt() << " hwQual "
0603                                     << result.muonCand->hwQual() << " hwEta " << result.muonCand->hwEta()
0604                                     << " deltaEta " << std::setw(8) << result.deltaEta << " deltaPhi " << std::setw(8)
0605                                     << result.deltaPhi << " Likelihood " << std::setw(8) << result.matchingLikelihood
0606                                     << " result " << (short)result.result;
0607       LogTrace("l1tOmtfEventPrint") << " procMuon " << *(result.procMuon) << std::endl;
0608     } else
0609       LogTrace("l1tOmtfEventPrint") << " no muonCand "
0610                                     << " result " << (short)result.result << std::endl;
0611   }
0612   LogTrace("l1tOmtfEventPrint") << " " << std::endl;
0613 
0614   return cleanedMatchingResults;
0615 }
0616 
0617 std::vector<MatchingResult> CandidateSimMuonMatcher::match(std::vector<const l1t::RegionalMuonCand*>& muonCands,
0618                                                            AlgoMuons& ghostBustedProcMuons,
0619                                                            const edm::SimTrackContainer* simTracks,
0620                                                            const edm::SimVertexContainer* simVertices,
0621                                                            std::function<bool(const SimTrack&)> const& simTrackFilter) {
0622   std::vector<MatchingResult> matchingResults;
0623 
0624   for (auto& simTrack : *simTracks) {
0625     if (!simTrackFilter(simTrack))
0626       continue;
0627 
0628     LogTrace("l1tOmtfEventPrint") << "CandidateSimMuonMatcher::match, simTrack type " << std::setw(3) << simTrack.type()
0629                                   << " pt " << std::setw(9) << simTrack.momentum().pt() << " eta " << std::setw(9)
0630                                   << simTrack.momentum().eta() << " phi " << std::setw(9) << simTrack.momentum().phi()
0631                                   << std::endl;
0632 
0633     bool matched = false;
0634 
0635     TrajectoryStateOnSurface tsof = propagate(simTrack, simVertices);
0636     if (!tsof.isValid()) {  //no sense to do matching
0637       MatchingResult result(simTrack);
0638       result.result = MatchingResult::ResultType::propagationFailed;
0639       LogTrace("l1tOmtfEventPrint") << __FUNCTION__ << ":" << __LINE__ << " propagation failed: genPt " << result.genPt
0640                                     << " genEta " << result.genEta << " eventId " << simTrack.eventId().event()
0641                                     << std::endl;
0642 
0643       matchingResults.push_back(result);
0644     } else {
0645       //checking if the propagated track is inside the OMTF range, TODO - tune the range!!!!!!!!!!!!!!!!!
0646       //eta 0.7 is the beginning of the MB2,
0647       //the eta range wider than the nominal OMTF region is needed, as in any case muons outside this region are seen by the OMTF
0648       //so it better to train the nn suich that is able to measure its pt, as it may affect the rate
0649       if ((std::abs(tsof.globalPosition().eta()) >= 0.7) && (std::abs(tsof.globalPosition().eta()) <= 1.31)) {
0650         LogTrace("l1tOmtfEventPrint")
0651             << "CandidateSimMuonMatcher::match simTrack IS in OMTF region, matching to the omtfCands";
0652       } else {
0653         LogTrace("l1tOmtfEventPrint") << "simTrack NOT in OMTF region ";
0654         continue;
0655       }
0656 
0657       /* TODO fix if filling of the deltaPhiPropCandMean and deltaPhiPropCandStdDev is needed
0658     double ptGen = simTrack.momentum().pt();
0659     if(ptGen >= deltaPhiVertexProp->GetXaxis()->GetXmax())
0660       ptGen = deltaPhiVertexProp->GetXaxis()->GetXmax() - 0.01;
0661 
0662     deltaPhiVertexProp->Fill(ptGen, simTrack.momentum().phi() - tsof.globalPosition().phi());*/
0663 
0664       unsigned int iCand = 0;
0665       for (auto& muonCand : muonCands) {
0666         //dropping very low quality candidates, as they are fakes usually - but it has no sense, then the results are not conclusive
0667         //if(muonCand->hwQual() > 1)
0668         {
0669           MatchingResult result;
0670           if (tsof.isValid()) {
0671             result = match(muonCand, ghostBustedProcMuons.at(iCand), simTrack, tsof);
0672           }
0673           int vtxInd = simTrack.vertIndex();
0674           if (vtxInd >= 0) {
0675             result.simVertex = &(simVertices->at(
0676                 vtxInd));  //TODO ?????? something strange is here, was commented in the previous version
0677           }
0678           if (result.result == MatchingResult::ResultType::matched) {
0679             matchingResults.push_back(result);
0680             matched = true;
0681           }
0682         }
0683         iCand++;
0684       }
0685 
0686       if (!matched) {  //we are adding also if it was not matching to any candidate
0687         MatchingResult result(simTrack);
0688         matchingResults.push_back(result);
0689         LogTrace("l1tOmtfEventPrint") << __FUNCTION__ << ":" << __LINE__ << " no matching candidate found" << std::endl;
0690       }
0691     }
0692   }
0693 
0694   return cleanMatching(matchingResults, muonCands, ghostBustedProcMuons);
0695 }
0696 
0697 std::vector<MatchingResult> CandidateSimMuonMatcher::match(
0698     std::vector<const l1t::RegionalMuonCand*>& muonCands,
0699     AlgoMuons& ghostBustedProcMuons,
0700     const TrackingParticleCollection* trackingParticles,
0701     std::function<bool(const TrackingParticle&)> const& simTrackFilter) {
0702   std::vector<MatchingResult> matchingResults;
0703   LogTrace("l1tOmtfEventPrint") << "CandidateSimMuonMatcher::match trackingParticles->size() "
0704                                 << trackingParticles->size() << std::endl;
0705 
0706   for (auto& trackingParticle : *trackingParticles) {
0707     //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;
0708 
0709     if (simTrackFilter(trackingParticle) == false)
0710       continue;
0711 
0712     LogTrace("l1tOmtfEventPrint") << "CandidateSimMuonMatcher::match, trackingParticle type " << std::setw(3)
0713                                   << trackingParticle.pdgId() << " pt " << std::setw(9) << trackingParticle.pt()
0714                                   << " eta " << std::setw(9) << trackingParticle.momentum().eta() << " phi "
0715                                   << std::setw(9) << trackingParticle.momentum().phi() << std::endl;
0716 
0717     bool matched = false;
0718 
0719     TrajectoryStateOnSurface tsof = propagate(trackingParticle);
0720     if (!tsof.isValid()) {
0721       LogTrace("l1tOmtfEventPrint") << "CandidateSimMuonMatcher::match:" << __LINE__ << " propagation failed"
0722                                     << std::endl;
0723       MatchingResult result;
0724       result.result = MatchingResult::ResultType::propagationFailed;
0725       continue;  //no sense to do matching
0726     }
0727 
0728     LogTrace("l1tOmtfEventPrint") << "CandidateSimMuonMatcher::match, tsof.globalPosition().eta() "
0729                                   << tsof.globalPosition().eta();
0730 
0731     //checking if the propagated track is inside the OMTF range, TODO - tune the range!!!!!!!!!!!!!!!!!
0732     //eta 0.7 is the beginning of the MB2, while 1.31 is mid of RE2 + some margin
0733     if ((std::abs(tsof.globalPosition().eta()) >= 0.7) && (std::abs(tsof.globalPosition().eta()) <= 1.31)) {
0734       LogTrace("l1tOmtfEventPrint")
0735           << "CandidateSimMuonMatcher::match trackingParticle IS in OMTF region, matching to the omtfCands";
0736     } else {
0737       LogTrace("l1tOmtfEventPrint") << "trackingParticle NOT in OMTF region ";
0738       continue;
0739     }
0740 
0741     /* TODO fix if filling of the deltaPhiPropCandMean and deltaPhiPropCandStdDev is needed
0742     double ptGen = trackingParticle.pt();
0743     if(ptGen >= deltaPhiVertexProp->GetXaxis()->GetXmax())
0744       ptGen = deltaPhiVertexProp->GetXaxis()->GetXmax() - 0.01;
0745 
0746     deltaPhiVertexProp->Fill(ptGen, trackingParticle.momentum().phi() - tsof.globalPosition().phi());
0747      */
0748 
0749     unsigned int iCand = 0;
0750     for (auto& muonCand : muonCands) {
0751       //dropping very low quality candidates, as they are fakes usually - but it has no sense, then the results are not conclusive then
0752       /*if(muonCand->hwQual() <= 1)
0753         continue; */
0754 
0755       MatchingResult result;
0756       if (tsof.isValid()) {
0757         result = match(muonCand, ghostBustedProcMuons.at(iCand), trackingParticle, tsof);
0758       }
0759       iCand++;
0760 
0761       if (result.result == MatchingResult::ResultType::matched) {
0762         matchingResults.push_back(result);
0763         matched = true;
0764       }
0765     }
0766 
0767     if (!matched) {  //we are adding result also if it there was no matching to any candidate
0768       MatchingResult result(trackingParticle);
0769       matchingResults.push_back(result);
0770       LogTrace("l1tOmtfEventPrint") << __FUNCTION__ << ":" << __LINE__ << " no matching candidate found" << std::endl;
0771     }
0772   }
0773 
0774   return cleanMatching(matchingResults, muonCands, ghostBustedProcMuons);
0775 }
0776 
0777 std::vector<MatchingResult> CandidateSimMuonMatcher::matchSimple(
0778     std::vector<const l1t::RegionalMuonCand*>& muonCands,
0779     AlgoMuons& ghostBustedProcMuons,
0780     const edm::SimTrackContainer* simTracks,
0781     const edm::SimVertexContainer* simVertices,
0782     std::function<bool(const SimTrack&)> const& simTrackFilter) {
0783   std::vector<MatchingResult> matchingResults;
0784 
0785   for (auto& simTrack : *simTracks) {
0786     if (!simTrackFilter(simTrack))
0787       continue;
0788 
0789     LogTrace("l1tOmtfEventPrint") << "CandidateSimMuonMatcher::matchSimple: simTrack type " << std::setw(3)
0790                                   << simTrack.type() << " pt " << std::setw(9) << simTrack.momentum().pt() << " eta "
0791                                   << std::setw(9) << simTrack.momentum().eta() << " phi " << std::setw(9)
0792                                   << simTrack.momentum().phi() << std::endl;
0793 
0794     bool matched = false;
0795 
0796     unsigned int iCand = 0;
0797     for (auto& muonCand : muonCands) {
0798       //dropping very low quality candidates, as they are fakes usually - but it has no sense, then the results are not conclusive
0799       //if(muonCand->hwQual() > 1)
0800       {
0801         MatchingResult result(simTrack);
0802 
0803         double candGloablEta = muonCand->hwEta() * 0.010875;
0804         double candGlobalPhi = omtfConfig->calcGlobalPhi(muonCand->hwPhi(), muonCand->processor());
0805         candGlobalPhi = hwGmtPhiToGlobalPhi(candGlobalPhi);
0806 
0807         if (candGlobalPhi > M_PI)
0808           candGlobalPhi = candGlobalPhi - (2. * M_PI);
0809 
0810         result.deltaPhi = foldPhi(result.genPhi - candGlobalPhi);
0811         result.deltaEta = result.genEta - candGloablEta;
0812 
0813         result.propagatedPhi = 0;
0814         result.propagatedEta = 0;
0815 
0816         result.muonCand = muonCand;
0817         result.procMuon = ghostBustedProcMuons.at(iCand);
0818 
0819         //TODO histogram can be used, like in the usercode/L1MuonAnalyzer  MuonMatcher::matchWithoutPorpagation
0820         //for prompt muons
0821         /*double treshold = 0.3;
0822         if (simTrack.momentum().pt() < 5)
0823           treshold = 1.5;
0824         else if (simTrack.momentum().pt() < 8)
0825           treshold = 0.8;
0826         else if (simTrack.momentum().pt() < 10)
0827           treshold = 0.6;
0828         else if (simTrack.momentum().pt() < 20)
0829           treshold = 0.5;*/
0830 
0831         //for displaced muons
0832         double treshold = 0.7;
0833         if (simTrack.momentum().pt() < 5)
0834           treshold = 1.5;
0835         else if (simTrack.momentum().pt() < 10)
0836           treshold = 1.0;
0837         else if (simTrack.momentum().pt() < 20)
0838           treshold = 0.7;
0839 
0840         if (std::abs(result.deltaPhi) < treshold && std::abs(result.deltaEta) < 0.5) {
0841           result.result = MatchingResult::ResultType::matched;
0842           //matchingLikelihood is needed in the cleanMatching, so we put something
0843           if (std::abs(result.deltaPhi) < 0.001)
0844             result.matchingLikelihood = 1. / 0.001;
0845           else
0846             result.matchingLikelihood = 1. / std::abs(result.deltaPhi);
0847         }
0848 
0849         LogTrace("l1tOmtfEventPrint") << "CandidateSimMuonMatcher::matchSimple: simTrack type " << simTrack.type()
0850                                       << " pt " << std::setw(8) << simTrack.momentum().pt() << " eta " << std::setw(8)
0851                                       << simTrack.momentum().eta() << " phi " << std::setw(8)
0852                                       << simTrack.momentum().phi() << "\n             muonCand pt " << std::setw(8)
0853                                       << muonCand->hwPt() << " candGloablEta " << std::setw(8) << candGloablEta
0854                                       << " candGlobalPhi " << std::setw(8) << candGlobalPhi << " hwQual "
0855                                       << muonCand->hwQual() << " deltaEta " << std::setw(8) << result.deltaEta
0856                                       << " deltaPhi " << std::setw(8) << result.deltaPhi << " matchingLikelihood "
0857                                       << result.matchingLikelihood << " result " << (short)result.result << std::endl;
0858 
0859         int vtxInd = simTrack.vertIndex();
0860         if (vtxInd >= 0) {
0861           result.simVertex = &(
0862               simVertices->at(vtxInd));  //TODO ?????? something strange is here, was commented in the previous version
0863         }
0864         if (result.result == MatchingResult::ResultType::matched) {
0865           matchingResults.push_back(result);
0866           matched = true;
0867         }
0868       }
0869       iCand++;
0870     }
0871 
0872     if (!matched) {  //we are adding also if it was not matched to any candidate
0873       MatchingResult result(simTrack);
0874       matchingResults.push_back(result);
0875       LogTrace("l1tOmtfEventPrint") << __FUNCTION__ << ":" << __LINE__ << " no matching candidate found" << std::endl;
0876     }
0877   }
0878 
0879   return cleanMatching(matchingResults, muonCands, ghostBustedProcMuons);
0880 }