File indexing completed on 2024-10-08 05:11:49
0001
0002
0003
0004
0005
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
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
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) {
0103
0104 } else
0105 return false;
0106
0107
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
0117
0118
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) {
0139
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) {
0150
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
0163
0164 } else
0165 return false;
0166
0167
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
0189
0190
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
0225
0226
0227 std::function<bool(const SimTrack&)> const& simTrackFilter = simTrackIsMuonInBx0;
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
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
0252 std::function<bool(const TrackingParticle&)> trackParticleFilter =
0253 trackingParticleIsMuonInBx0;
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
0288 int deltaPhi = std::abs(gloablHwPhi1 - gloablHwPhi2);
0289 if (deltaPhi > 576 / 2)
0290 deltaPhi = std::abs(deltaPhi - 576);
0291
0292
0293
0294
0295
0296
0297 if (deltaPhi < 8) {
0298
0299
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
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
0344
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;
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;
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);
0391
0392 return tsof;
0393 }
0394
0395 TrajectoryStateOnSurface CandidateSimMuonMatcher::propagate(const TrackingParticle& trackingParticle) {
0396 FreeTrajectoryState ftsTrack = simTrackToFts(trackingParticle);
0397
0398 TrajectoryStateOnSurface tsof = atStation2(ftsTrack);
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
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
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);
0440
0441 result.muonCand = muonCand;
0442 result.procMuon = procMuon;
0443
0444
0445 double treshold = 0.15;
0446 if (simTrack.momentum().pt() <
0447 10)
0448 treshold = 0.3;
0449 else if (simTrack.momentum().pt() <
0450 30)
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
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
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);
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
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
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
0560 cleanedMatchingResults.push_back(matchingResult);
0561 if (matchingResult.result == MatchingResult::ResultType::matched) {
0562
0563
0564
0565
0566
0567
0568 }
0569 }
0570
0571
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()) {
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
0646
0647
0648
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
0658
0659
0660
0661
0662
0663
0664 unsigned int iCand = 0;
0665 for (auto& muonCand : muonCands) {
0666
0667
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));
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) {
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
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;
0726 }
0727
0728 LogTrace("l1tOmtfEventPrint") << "CandidateSimMuonMatcher::match, tsof.globalPosition().eta() "
0729 << tsof.globalPosition().eta();
0730
0731
0732
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
0742
0743
0744
0745
0746
0747
0748
0749 unsigned int iCand = 0;
0750 for (auto& muonCand : muonCands) {
0751
0752
0753
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) {
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
0799
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
0820
0821
0822
0823
0824
0825
0826
0827
0828
0829
0830
0831
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
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));
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) {
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 }