Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 11:12:54 

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 "TFile.h"
0024 #include "TH1D.h"
0025 
0026 //#include "CLHEP/Units/GlobalPhysicalConstants.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   deltaPhiPropCandMean = (TH1D*)inFile.Get("deltaPhiPropCandMean");
0057   deltaPhiPropCandStdDev = (TH1D*)inFile.Get("deltaPhiPropCandStdDev");
0058 }
0059 
0060 CandidateSimMuonMatcher::~CandidateSimMuonMatcher() {}
0061 
0062 void CandidateSimMuonMatcher::beginRun(const edm::EventSetup& eventSetup) {
0063   //TODO use edm::ESWatcher<MagneticField> magneticFieldRecordWatcher;
0064   magField = eventSetup.getHandle(magneticFieldEsToken);
0065   propagator = eventSetup.getHandle(propagatorEsToken);
0066 }
0067 
0068 void CandidateSimMuonMatcher::observeEventBegin(const edm::Event& event) { gbCandidates.clear(); }
0069 
0070 void CandidateSimMuonMatcher::observeProcesorEmulation(unsigned int iProcessor,
0071                                                        l1t::tftype mtfType,
0072                                                        const std::shared_ptr<OMTFinput>& input,
0073                                                        const AlgoMuons& algoCandidates,
0074                                                        const AlgoMuons& gbCandidates,
0075                                                        const std::vector<l1t::RegionalMuonCand>& candMuons) {
0076   //debug
0077   /*
0078   unsigned int procIndx = omtfConfig->getProcIndx(iProcessor, mtfType);
0079   for (auto& gbCandidate : gbCandidates) {
0080     if (gbCandidate->getPt() > 0) {
0081       LogTrace("l1tOmtfEventPrint") << "CandidateSimMuonMatcher::observeProcesorEmulation procIndx" << procIndx<" "<< *gbCandidate << endl;
0082       this->gbCandidates.emplace_back(gbCandidate);
0083     }
0084   }*/
0085 }
0086 
0087 bool simTrackIsMuonInOmtf(const SimTrack& simTrack) {
0088   if (abs(simTrack.type()) == 13 || abs(simTrack.type()) == 1000015) {  //|| tpPtr->pt() > 20 //todo 1000015 is stau
0089     //only muons
0090   } else
0091     return false;
0092 
0093   //in the overlap, the propagation of muons with pt less then ~3.2 fails - the actual threshold depends slightly on eta,
0094   if (simTrack.momentum().pt() < 2.5)
0095     return false;
0096 
0097   LogTrace("l1tOmtfEventPrint") << "simTrackIsMuonInOmtf, simTrack type " << std::setw(3) << simTrack.type() << " pt "
0098                                 << std::setw(9) << simTrack.momentum().pt() << " eta " << std::setw(9)
0099                                 << simTrack.momentum().eta() << " phi " << std::setw(9) << simTrack.momentum().phi()
0100                                 << std::endl;
0101 
0102   //higher margin for matching must be used than actual OMTF region (i.e. 0.82 - 1.24),
0103   //otherwise many candidates are marked as ghosts
0104   //if( (fabs(simTrack.momentum().eta()) >= 0.82 ) && (fabs(simTrack.momentum().eta()) <= 1.24) ) {
0105   if ((fabs(simTrack.momentum().eta()) >= 0.72) && (fabs(simTrack.momentum().eta()) <= 1.3)) {
0106   } else
0107     return false;
0108 
0109   return true;
0110 }
0111 
0112 bool simTrackIsMuonInOmtfBx0(const SimTrack& simTrack) {
0113   if (simTrack.eventId().bunchCrossing() != 0)
0114     return false;
0115 
0116   return simTrackIsMuonInOmtf(simTrack);
0117 }
0118 
0119 bool trackingParticleIsMuonInOmtfBx0(const TrackingParticle& trackingParticle) {
0120   /*if(trackingParticle.eventId().event() != 0)
0121       LogTrace("l1tOmtfEventPrint") <<"trackingParticleIsMuonInOmtfBx0, pdgId "<<std::setw(3)<<trackingParticle.pdgId()<<" pt "<<std::setw(9)<<trackingParticle.pt()
0122               <<" eta "<<std::setw(9)<<trackingParticle.momentum().eta()<<" phi "<<std::setw(9)<<trackingParticle.momentum().phi()<<" event "<<trackingParticle.eventId().event()
0123               <<" bx "<<trackingParticle.eventId().bunchCrossing()<<" eventNot0"<<std::endl;*/
0124 
0125   if (trackingParticle.eventId().bunchCrossing() != 0)
0126     return false;
0127 
0128   if (abs(trackingParticle.pdgId()) == 13 || abs(trackingParticle.pdgId()) == 1000015) {
0129     // 1000015 is stau, todo use other selection if needed
0130     //only muons
0131   } else
0132     return false;
0133 
0134   //in the overlap, the propagation of muons with pt less then ~3.2 fails, the actual threshold depends slightly on eta,
0135   if (trackingParticle.pt() < 2.5)
0136     return false;
0137 
0138   if (trackingParticle.parentVertex().isNonnull())
0139     LogTrace("l1tOmtfEventPrint") << "trackingParticleIsMuonInOmtfBx0, pdgId " << std::setw(3)
0140                                   << trackingParticle.pdgId() << " pt " << std::setw(9) << trackingParticle.pt()
0141                                   << " eta " << std::setw(9) << trackingParticle.momentum().eta() << " phi "
0142                                   << std::setw(9) << trackingParticle.momentum().phi() << " event "
0143                                   << trackingParticle.eventId().event() << " trackId "
0144                                   << trackingParticle.g4Tracks().at(0).trackId() << " parentVertex Rho "
0145                                   << trackingParticle.parentVertex()->position().Rho() << " eta "
0146                                   << trackingParticle.parentVertex()->position().eta() << " phi "
0147                                   << trackingParticle.parentVertex()->position().phi() << std::endl;
0148   else
0149     LogTrace("l1tOmtfEventPrint") << "trackingParticleIsMuonInOmtfBx0, pdgId " << std::setw(3)
0150                                   << trackingParticle.pdgId() << " pt " << std::setw(9) << trackingParticle.pt()
0151                                   << " eta " << std::setw(9) << trackingParticle.momentum().eta() << " phi "
0152                                   << std::setw(9) << trackingParticle.momentum().phi() << " trackId "
0153                                   << trackingParticle.g4Tracks().at(0).trackId();
0154 
0155   //higher margin for matching must be used than actual OMTF region (i.e. 0.82 - 1.24),
0156   //otherwise many candidates are marked as ghosts
0157   //if( (fabs(simTrack.momentum().eta()) >= 0.82 ) && (fabs(trackingParticle.momentum().eta()) <= 1.24) ) {
0158   if ((fabs(trackingParticle.momentum().eta()) >= 0.72) && (fabs(trackingParticle.momentum().eta()) <= 1.3)) {
0159   } else
0160     return false;
0161 
0162   return true;
0163 }
0164 
0165 bool trackingParticleIsMuonInOmtfEvent0(const TrackingParticle& trackingParticle) {
0166   if (trackingParticle.eventId().event() != 0)
0167     return false;
0168 
0169   return trackingParticleIsMuonInOmtfBx0(trackingParticle);
0170 }
0171 
0172 void CandidateSimMuonMatcher::observeEventEnd(const edm::Event& event,
0173                                               std::unique_ptr<l1t::RegionalMuonCandBxCollection>& finalCandidates) {
0174   AlgoMuons ghostBustedProcMuons;
0175   std::vector<const l1t::RegionalMuonCand*> ghostBustedRegionalCands =
0176       CandidateSimMuonMatcher::ghostBust(finalCandidates.get(), gbCandidates, ghostBustedProcMuons);
0177 
0178   matchingResults.clear();
0179   if (edmCfg.exists("simTracksTag")) {
0180     edm::Handle<edm::SimTrackContainer> simTraksHandle;
0181     event.getByLabel(edmCfg.getParameter<edm::InputTag>("simTracksTag"), simTraksHandle);
0182 
0183     edm::Handle<edm::SimVertexContainer> simVertices;
0184     event.getByLabel(edmCfg.getParameter<edm::InputTag>("simVertexesTag"), simVertices);
0185     LogTrace("l1tOmtfEventPrint") << "simTraksHandle size " << simTraksHandle.product()->size() << std::endl;
0186 
0187     //TODO  use other simTrackFilter if needed  <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
0188     std::function<bool(const SimTrack&)> const& simTrackFilter = simTrackIsMuonInOmtfBx0;
0189 
0190     matchingResults = match(
0191         ghostBustedRegionalCands, ghostBustedProcMuons, simTraksHandle.product(), simVertices.product(), simTrackFilter);
0192 
0193   } else if (edmCfg.exists("trackingParticleTag")) {
0194     edm::Handle<TrackingParticleCollection> trackingParticleHandle;
0195     event.getByLabel(edmCfg.getParameter<edm::InputTag>("trackingParticleTag"), trackingParticleHandle);
0196     LogTrace("l1tOmtfEventPrint") << "trackingParticleHandle size " << trackingParticleHandle.product()->size()
0197                                   << std::endl;
0198 
0199     //TODO use other trackParticleFilter if needed <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
0200     std::function<bool(const TrackingParticle&)> trackParticleFilter = trackingParticleIsMuonInOmtfBx0;
0201     matchingResults =
0202         match(ghostBustedRegionalCands, ghostBustedProcMuons, trackingParticleHandle.product(), trackParticleFilter);
0203   }
0204 }
0205 
0206 void CandidateSimMuonMatcher::endJob() {}
0207 
0208 std::vector<const l1t::RegionalMuonCand*> CandidateSimMuonMatcher::ghostBust(
0209     const l1t::RegionalMuonCandBxCollection* mtfCands, const AlgoMuons& gbCandidates, AlgoMuons& ghostBustedProcMuons) {
0210   if (gbCandidates.size() != mtfCands->size()) {
0211     edm::LogError("l1tOmtfEventPrint") << "CandidateSimMuonMatcher::ghostBust(): gbCandidates.size() "
0212                                        << gbCandidates.size() << " != resultGbCandidates.size() " << mtfCands->size();
0213     throw cms::Exception("gbCandidates.size() != mtfCands->size()");
0214   }
0215 
0216   boost::dynamic_bitset<> isKilled(mtfCands->size(0), false);
0217 
0218   for (unsigned int i1 = 0; i1 < mtfCands->size(0); ++i1) {
0219     for (unsigned int i2 = i1 + 1; i2 < mtfCands->size(0); ++i2) {
0220       auto& mtfCand1 = mtfCands->at(0, i1);
0221       auto& mtfCand2 = mtfCands->at(0, i2);
0222 
0223       if (abs(mtfCand1.hwEta() - mtfCand2.hwEta()) < (0.3 / 0.010875)) {
0224         int gloablHwPhi1 = l1t::MicroGMTConfiguration::calcGlobalPhi(
0225             mtfCand1.hwPhi(), mtfCand1.trackFinderType(), mtfCand1.processor());
0226         int gloablHwPhi2 = l1t::MicroGMTConfiguration::calcGlobalPhi(
0227             mtfCand2.hwPhi(), mtfCand2.trackFinderType(), mtfCand2.processor());
0228 
0229         //one can use the phi in radians like that:
0230         //double globalPhi1 = hwGmtPhiToGlobalPhi(l1t::MicroGMTConfiguration::calcGlobalPhi( mtfCand1.hwPhi(), mtfCand1.trackFinderType(), mtfCand1.processor() ) );
0231         //double globalPhi2 = hwGmtPhiToGlobalPhi(l1t::MicroGMTConfiguration::calcGlobalPhi( mtfCand2.hwPhi(), mtfCand2.trackFinderType(), mtfCand2.processor() ) );
0232 
0233         //0.0872664626 = 5 deg, i.e. the same window as in the OMTF ghost buster
0234         if (abs(gloablHwPhi1 - gloablHwPhi2) < 8) {
0235           if (mtfCand1.hwQual() > mtfCand2.hwQual()) {
0236             isKilled[i2] = true;
0237           } else
0238             isKilled[i1] = true;
0239         }
0240       }
0241     }
0242   }
0243 
0244   std::vector<const l1t::RegionalMuonCand*> resultCands;
0245 
0246   for (unsigned int i1 = 0; i1 < mtfCands->size(0); ++i1) {
0247     //dropping candidates with quality 0 !!!!!!!!!!!!!!!!!!!! fixme if not needed
0248     if (!isKilled[i1] && mtfCands->at(0, i1).hwQual()) {
0249       resultCands.push_back(&(mtfCands->at(0, i1)));
0250       ghostBustedProcMuons.push_back(gbCandidates.at(i1));
0251     }
0252 
0253     LogTrace("l1tOmtfEventPrint")
0254         << "CandidateSimMuonMatcher::ghostBust\n  regionalCand pt " << std::setw(3) << mtfCands->at(0, i1).hwPt()
0255         << " qual " << std::setw(2) << mtfCands->at(0, i1).hwQual() << " proc " << std::setw(2)
0256         << mtfCands->at(0, i1).processor() << " eta " << std::setw(4) << mtfCands->at(0, i1).hwEta() << " gloablEta "
0257         << std::setw(8) << mtfCands->at(0, i1).hwEta() * 0.010875 << " hwPhi " << std::setw(3)
0258         << mtfCands->at(0, i1).hwPhi() << " globalPhi " << std::setw(8)
0259         << hwGmtPhiToGlobalPhi(l1t::MicroGMTConfiguration::calcGlobalPhi(
0260                mtfCands->at(0, i1).hwPhi(), mtfCands->at(0, i1).trackFinderType(), mtfCands->at(0, i1).processor()))
0261         << " fireadLayers " << std::bitset<18>(mtfCands->at(0, i1).trackAddress().at(0)) << " isKilled "
0262         << isKilled.test(i1) << std::endl;
0263 
0264     LogTrace("l1tOmtfEventPrint") << *(gbCandidates.at(i1)) << std::endl;
0265   }
0266 
0267   if (resultCands.size() >= 3)
0268     LogTrace("l1tOmtfEventPrint") << " ghost !!!!!! " << std::endl;
0269   LogTrace("l1tOmtfEventPrint") << std::endl;
0270 
0271   return resultCands;
0272 }
0273 
0274 TrajectoryStateOnSurface CandidateSimMuonMatcher::atStation2(FreeTrajectoryState ftsStart, float eta) const {
0275   ReferenceCountingPointer<Surface> rpc;
0276   if (eta < -1.24)  //negative endcap, RE2
0277     rpc = ReferenceCountingPointer<Surface>(
0278         new BoundDisk(GlobalPoint(0., 0., -790.), TkRotation<float>(), SimpleDiskBounds(300., 810., -10., 10.)));
0279   else if (eta < 1.24)  //barrel + overlap, 512.401cm is R of middle of the MB2
0280     rpc = ReferenceCountingPointer<Surface>(new BoundCylinder(
0281         GlobalPoint(0., 0., 0.), TkRotation<float>(), SimpleCylinderBounds(512.401, 512.401, -900, 900)));
0282   else
0283     rpc = ReferenceCountingPointer<Surface>(  //positive endcap, RE2
0284         new BoundDisk(GlobalPoint(0., 0., 790.), TkRotation<float>(), SimpleDiskBounds(300., 810., -10., 10.)));
0285 
0286   TrajectoryStateOnSurface trackAtRPC = propagator->propagate(ftsStart, *rpc);
0287   return trackAtRPC;
0288 }
0289 
0290 FreeTrajectoryState CandidateSimMuonMatcher::simTrackToFts(const SimTrack& simTrackPtr, const SimVertex& simVertex) {
0291   int charge = simTrackPtr.type() > 0 ? -1 : 1;  //works for muons
0292 
0293   GlobalVector p3GV(simTrackPtr.momentum().x(), simTrackPtr.momentum().y(), simTrackPtr.momentum().z());
0294   GlobalPoint r3GP(simVertex.position().x(), simVertex.position().y(), simVertex.position().z());
0295 
0296   GlobalTrajectoryParameters tPars(r3GP, p3GV, charge, &*magField);
0297 
0298   return FreeTrajectoryState(tPars);
0299 }
0300 
0301 FreeTrajectoryState CandidateSimMuonMatcher::simTrackToFts(const TrackingParticle& trackingParticle) {
0302   int charge = trackingParticle.pdgId() > 0 ? -1 : 1;  //works for muons
0303 
0304   GlobalVector p3GV(trackingParticle.momentum().x(), trackingParticle.momentum().y(), trackingParticle.momentum().z());
0305   GlobalPoint r3GP(trackingParticle.vx(), trackingParticle.vy(), trackingParticle.vz());
0306 
0307   GlobalTrajectoryParameters tPars(r3GP, p3GV, charge, &*magField);
0308 
0309   return FreeTrajectoryState(tPars);
0310 }
0311 
0312 TrajectoryStateOnSurface CandidateSimMuonMatcher::propagate(const SimTrack& simTrack,
0313                                                             const edm::SimVertexContainer* simVertices) {
0314   SimVertex simVertex;
0315   int vtxInd = simTrack.vertIndex();
0316   if (vtxInd < 0) {
0317     std::cout << "Track with no vertex, defaulting to (0,0,0)" << std::endl;
0318   } else {
0319     simVertex = simVertices->at(vtxInd);
0320     if (((int)simVertex.vertexId()) != vtxInd) {
0321       std::cout << "simVertex.vertexId() != vtxInd !!!!!!!!!!!!!!!!!" << std::endl;
0322       edm::LogImportant("l1tOmtfEventPrint") << "simVertex.vertexId() != vtxInd. simVertex.vertexId() "
0323                                              << simVertex.vertexId() << " vtxInd " << vtxInd << " !!!!!!!!!!!!!!!!!";
0324     }
0325   }
0326 
0327   FreeTrajectoryState ftsTrack = simTrackToFts(simTrack, simVertex);
0328 
0329   TrajectoryStateOnSurface tsof = atStation2(ftsTrack, simTrack.momentum().eta());  //propagation
0330 
0331   return tsof;
0332 }
0333 
0334 TrajectoryStateOnSurface CandidateSimMuonMatcher::propagate(const TrackingParticle& trackingParticle) {
0335   FreeTrajectoryState ftsTrack = simTrackToFts(trackingParticle);
0336 
0337   TrajectoryStateOnSurface tsof = atStation2(ftsTrack, trackingParticle.momentum().eta());  //propagation
0338 
0339   return tsof;
0340 }
0341 
0342 float normal_pdf(float x, float m, float s) {
0343   static const float inv_sqrt_2pi = 0.3989422804014327;
0344   float a = (x - m) / s;
0345 
0346   return inv_sqrt_2pi / s * std::exp(-0.5 * a * a);
0347 }
0348 
0349 MatchingResult CandidateSimMuonMatcher::match(const l1t::RegionalMuonCand* muonCand,
0350                                               const AlgoMuonPtr& procMuon,
0351                                               const SimTrack& simTrack,
0352                                               TrajectoryStateOnSurface& tsof) {
0353   MatchingResult result(simTrack);
0354 
0355   double candGloablEta = muonCand->hwEta() * 0.010875;
0356   if (abs(simTrack.momentum().eta() - candGloablEta) < 0.3) {
0357     double candGlobalPhi = l1t::MicroGMTConfiguration::calcGlobalPhi(
0358         muonCand->hwPhi(), muonCand->trackFinderType(), muonCand->processor());
0359     candGlobalPhi = hwGmtPhiToGlobalPhi(candGlobalPhi);
0360 
0361     if (candGlobalPhi > M_PI)
0362       candGlobalPhi = candGlobalPhi - (2. * M_PI);
0363 
0364     result.deltaPhi = foldPhi(tsof.globalPosition().phi() - candGlobalPhi);
0365     result.deltaEta = tsof.globalPosition().eta() - candGloablEta;
0366 
0367     result.propagatedPhi = tsof.globalPosition().phi();
0368     result.propagatedEta = tsof.globalPosition().eta();
0369 
0370     double mean = 0;
0371     double sigma = 1;
0372     //if(!fillMean)
0373     {
0374       auto ptBin = deltaPhiPropCandMean->FindBin(simTrack.momentum().pt());
0375       mean = deltaPhiPropCandMean->GetBinContent(ptBin);
0376       sigma = deltaPhiPropCandStdDev->GetBinContent(ptBin);
0377     }
0378     result.matchingLikelihood = normal_pdf(result.deltaPhi, mean, sigma);  //TODO temporary solution
0379 
0380     result.muonCand = muonCand;
0381     result.procMuon = procMuon;
0382 
0383     double treshold = 6. * sigma;
0384     if (simTrack.momentum().pt() > 20)
0385       treshold = 7. * sigma;
0386     if (simTrack.momentum().pt() > 100)
0387       treshold = 20. * sigma;
0388 
0389     if (abs(result.deltaPhi - mean) < treshold)
0390       result.result = MatchingResult::ResultType::matched;
0391 
0392     LogTrace("l1tOmtfEventPrint") << "CandidateSimMuonMatcher::match: simTrack type " << simTrack.type() << " pt "
0393                                   << std::setw(8) << simTrack.momentum().pt() << " eta " << std::setw(8)
0394                                   << simTrack.momentum().eta() << " phi " << std::setw(8) << simTrack.momentum().phi()
0395                                   << " propagation eta " << std::setw(8) << tsof.globalPosition().eta() << " phi "
0396                                   << tsof.globalPosition().phi() << "\n             muonCand pt " << std::setw(8)
0397                                   << muonCand->hwPt() << " candGloablEta " << std::setw(8) << candGloablEta
0398                                   << " candGlobalPhi " << std::setw(8) << candGlobalPhi << " hwQual "
0399                                   << muonCand->hwQual() << " deltaEta " << std::setw(8) << result.deltaEta
0400                                   << " deltaPhi " << std::setw(8) << result.deltaPhi << " sigma " << std::setw(8)
0401                                   << sigma << " Likelihood " << std::setw(8) << result.matchingLikelihood << " result "
0402                                   << (short)result.result << std::endl;
0403   }
0404 
0405   return result;
0406 }
0407 
0408 MatchingResult CandidateSimMuonMatcher::match(const l1t::RegionalMuonCand* muonCand,
0409                                               const AlgoMuonPtr& procMuon,
0410                                               const TrackingParticle& trackingParticle,
0411                                               TrajectoryStateOnSurface& tsof) {
0412   MatchingResult result(trackingParticle);
0413 
0414   double candGloablEta = muonCand->hwEta() * 0.010875;
0415   if (abs(trackingParticle.momentum().eta() - candGloablEta) < 0.3) {
0416     double candGlobalPhi = l1t::MicroGMTConfiguration::calcGlobalPhi(
0417         muonCand->hwPhi(), muonCand->trackFinderType(), muonCand->processor());
0418     candGlobalPhi = hwGmtPhiToGlobalPhi(candGlobalPhi);
0419 
0420     if (candGlobalPhi > M_PI)
0421       candGlobalPhi = candGlobalPhi - (2. * M_PI);
0422 
0423     result.deltaPhi = foldPhi(tsof.globalPosition().phi() - candGlobalPhi);
0424     result.deltaEta = tsof.globalPosition().eta() - candGloablEta;
0425 
0426     result.propagatedPhi = tsof.globalPosition().phi();
0427     result.propagatedEta = tsof.globalPosition().eta();
0428 
0429     double mean = 0;
0430     double sigma = 1;
0431     //if(!fillMean)
0432     {
0433       auto ptBin = deltaPhiPropCandMean->FindBin(trackingParticle.pt());
0434 
0435       mean = deltaPhiPropCandMean->GetBinContent(ptBin);
0436       sigma = deltaPhiPropCandStdDev->GetBinContent(ptBin);
0437     }
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 (trackingParticle.pt() > 20)
0446       treshold = 7. * sigma;
0447     if (trackingParticle.pt() > 100)
0448       treshold = 20. * sigma;
0449 
0450     if (abs(result.deltaPhi - mean) < treshold)
0451       result.result = MatchingResult::ResultType::matched;
0452 
0453     LogTrace("l1tOmtfEventPrint") << "CandidateSimMuonMatcher::match: simTrack type " << trackingParticle.pdgId()
0454                                   << " pt " << std::setw(8) << trackingParticle.pt() << " eta " << std::setw(8)
0455                                   << trackingParticle.momentum().eta() << " phi " << std::setw(8)
0456                                   << trackingParticle.momentum().phi() << " propagation eta " << std::setw(8)
0457                                   << tsof.globalPosition().eta() << " phi " << tsof.globalPosition().phi()
0458                                   << " muonCand pt " << std::setw(8) << muonCand->hwPt() << " candGloablEta "
0459                                   << std::setw(8) << candGloablEta << " candGlobalPhi " << std::setw(8) << candGlobalPhi
0460                                   << " hwQual " << muonCand->hwQual() << " deltaEta " << std::setw(8) << result.deltaEta
0461                                   << " deltaPhi " << std::setw(8) << result.deltaPhi << " Likelihood " << std::setw(8)
0462                                   << result.matchingLikelihood << " result " << (short)result.result << std::endl;
0463   }
0464 
0465   return result;
0466 }
0467 
0468 std::vector<MatchingResult> CandidateSimMuonMatcher::cleanMatching(std::vector<MatchingResult> matchingResults,
0469                                                                    std::vector<const l1t::RegionalMuonCand*>& muonCands,
0470                                                                    AlgoMuons& ghostBustedProcMuons) {
0471   //Cleaning the matching
0472   std::sort(
0473       matchingResults.begin(), matchingResults.end(), [](const MatchingResult& a, const MatchingResult& b) -> bool {
0474         return a.matchingLikelihood > b.matchingLikelihood;
0475       });
0476   for (unsigned int i1 = 0; i1 < matchingResults.size(); i1++) {
0477     if (matchingResults[i1].result == MatchingResult::ResultType::matched) {
0478       for (unsigned int i2 = i1 + 1; i2 < matchingResults.size(); i2++) {
0479         if ((matchingResults[i1].trackingParticle &&
0480              matchingResults[i1].trackingParticle == matchingResults[i2].trackingParticle) ||
0481             (matchingResults[i1].simTrack && matchingResults[i1].simTrack == matchingResults[i2].simTrack) ||
0482             (matchingResults[i1].muonCand == matchingResults[i2].muonCand)) {
0483           //if matchingResults[i1].muonCand == false, then it is also OK here
0484           matchingResults[i2].result = MatchingResult::ResultType::duplicate;
0485         }
0486       }
0487     }
0488   }
0489 
0490   std::vector<MatchingResult> cleanedMatchingResults;
0491   for (auto& matchingResult : matchingResults) {
0492     if (matchingResult.result == MatchingResult::ResultType::matched ||
0493         matchingResult.muonCand ==
0494             nullptr)  //adding also the simTracks that are not matched at all, before it is assured that they are not duplicates
0495       cleanedMatchingResults.push_back(matchingResult);
0496     if (matchingResult.result == MatchingResult::ResultType::matched) {
0497       /* TODO fix if filling of the deltaPhiPropCandMean and deltaPhiPropCandStdDev is needed
0498       if(fillMean) {
0499       double ptGen  = matchingResult.genPt;
0500         deltaPhiPropCandMean->Fill(ptGen, matchingResult.deltaPhi); //filling overflow is ok here
0501         deltaPhiPropCandStdDev->Fill(ptGen, matchingResult.deltaPhi * matchingResult.deltaPhi);
0502       }*/
0503     }
0504   }
0505 
0506   //adding the muonCand-s that were not matched, i.e. in order to analyze them later
0507   unsigned int iCand = 0;
0508   for (auto& muonCand : muonCands) {
0509     bool isMatched = false;
0510     for (auto& matchingResult : cleanedMatchingResults) {
0511       if (matchingResult.muonCand == muonCand) {
0512         isMatched = true;
0513         break;
0514       }
0515     }
0516 
0517     if (!isMatched) {
0518       MatchingResult result;
0519       result.muonCand = muonCand;
0520       result.procMuon = ghostBustedProcMuons.at(iCand);
0521       cleanedMatchingResults.push_back(result);
0522     }
0523     iCand++;
0524   }
0525 
0526   LogTrace("l1tOmtfEventPrint") << "CandidateSimMuonMatcher::cleanMatching:" << __LINE__
0527                                 << " CandidateSimMuonMatcher::match cleanedMatchingResults:" << std::endl;
0528   for (auto& result : cleanedMatchingResults) {
0529     if (result.trackingParticle || result.simTrack)
0530       LogTrace("l1tOmtfEventPrint") << "CandidateSimMuonMatcher::cleanMatching:" << __LINE__ << " simTrack type "
0531                                     << result.pdgId << " pt " << std::setw(8) << result.genPt << " eta " << std::setw(8)
0532                                     << result.genEta << " phi " << std::setw(8) << result.genPhi;
0533     else
0534       LogTrace("l1tOmtfEventPrint") << "no matched track ";
0535 
0536     if (result.muonCand) {
0537       LogTrace("l1tOmtfEventPrint") << " muonCand pt " << std::setw(8) << result.muonCand->hwPt() << " hwQual "
0538                                     << result.muonCand->hwQual() << " hwEta " << result.muonCand->hwEta()
0539                                     << " deltaEta " << std::setw(8) << result.deltaEta << " deltaPhi " << std::setw(8)
0540                                     << result.deltaPhi << " Likelihood " << std::setw(8) << result.matchingLikelihood
0541                                     << " result " << (short)result.result;
0542       LogTrace("l1tOmtfEventPrint") << " procMuon " << *(result.procMuon) << std::endl;
0543     } else
0544       LogTrace("l1tOmtfEventPrint") << " no muonCand "
0545                                     << " result " << (short)result.result << std::endl;
0546   }
0547   LogTrace("l1tOmtfEventPrint") << " " << std::endl;
0548 
0549   return cleanedMatchingResults;
0550 }
0551 
0552 std::vector<MatchingResult> CandidateSimMuonMatcher::match(std::vector<const l1t::RegionalMuonCand*>& muonCands,
0553                                                            AlgoMuons& ghostBustedProcMuons,
0554                                                            const edm::SimTrackContainer* simTracks,
0555                                                            const edm::SimVertexContainer* simVertices,
0556                                                            std::function<bool(const SimTrack&)> const& simTrackFilter) {
0557   std::vector<MatchingResult> matchingResults;
0558 
0559   for (auto& simTrack : *simTracks) {
0560     if (!simTrackFilter(simTrack))
0561       continue;
0562 
0563     LogTrace("l1tOmtfEventPrint") << "CandidateSimMuonMatcher::match, simTrack type " << std::setw(3) << simTrack.type()
0564                                   << " pt " << std::setw(9) << simTrack.momentum().pt() << " eta " << std::setw(9)
0565                                   << simTrack.momentum().eta() << " phi " << std::setw(9) << simTrack.momentum().phi()
0566                                   << std::endl;
0567 
0568     bool matched = false;
0569 
0570     TrajectoryStateOnSurface tsof = propagate(simTrack, simVertices);
0571     if (!tsof.isValid()) {
0572       LogTrace("l1tOmtfEventPrint") << __FUNCTION__ << ":" << __LINE__ << " propagation failed" << std::endl;
0573       MatchingResult result;
0574       result.result = MatchingResult::ResultType::propagationFailed;
0575       continue;  //no sense to do matching
0576     }
0577 
0578     /* TODO fix if filling of the deltaPhiPropCandMean and deltaPhiPropCandStdDev is needed
0579     double ptGen = simTrack.momentum().pt();
0580     if(ptGen >= deltaPhiVertexProp->GetXaxis()->GetXmax())
0581       ptGen = deltaPhiVertexProp->GetXaxis()->GetXmax() - 0.01;
0582 
0583     deltaPhiVertexProp->Fill(ptGen, simTrack.momentum().phi() - tsof.globalPosition().phi());*/
0584 
0585     unsigned int iCand = 0;
0586     for (auto& muonCand : muonCands) {
0587       //dropping very low quality candidates, as they are fakes usually - but it has no sense, then the results are not conclusive
0588       //if(muonCand->hwQual() > 1)
0589       {
0590         MatchingResult result = match(muonCand, ghostBustedProcMuons.at(iCand), simTrack, tsof);
0591         if (result.result == MatchingResult::ResultType::matched) {
0592           matchingResults.push_back(result);
0593           matched = true;
0594         }
0595       }
0596       iCand++;
0597     }
0598 
0599     if (!matched) {  //we are adding also if it was not matching to any candidate
0600       MatchingResult result(simTrack);
0601       matchingResults.push_back(result);
0602       LogTrace("l1tOmtfEventPrint") << __FUNCTION__ << ":" << __LINE__ << " no matching candidate found" << std::endl;
0603     }
0604   }
0605 
0606   return cleanMatching(matchingResults, muonCands, ghostBustedProcMuons);
0607 }
0608 
0609 std::vector<MatchingResult> CandidateSimMuonMatcher::match(
0610     std::vector<const l1t::RegionalMuonCand*>& muonCands,
0611     AlgoMuons& ghostBustedProcMuons,
0612     const TrackingParticleCollection* trackingParticles,
0613     std::function<bool(const TrackingParticle&)> const& simTrackFilter) {
0614   std::vector<MatchingResult> matchingResults;
0615   LogTrace("l1tOmtfEventPrint") << "CandidateSimMuonMatcher::match trackingParticles->size() "
0616                                 << trackingParticles->size() << std::endl;
0617 
0618   for (auto& trackingParticle : *trackingParticles) {
0619     //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;
0620 
0621     if (simTrackFilter(trackingParticle) == false)
0622       continue;
0623 
0624     LogTrace("l1tOmtfEventPrint") << "CandidateSimMuonMatcher::match, trackingParticle type " << std::setw(3)
0625                                   << trackingParticle.pdgId() << " pt " << std::setw(9) << trackingParticle.pt()
0626                                   << " eta " << std::setw(9) << trackingParticle.momentum().eta() << " phi "
0627                                   << std::setw(9) << trackingParticle.momentum().phi() << std::endl;
0628 
0629     bool matched = false;
0630 
0631     TrajectoryStateOnSurface tsof = propagate(trackingParticle);
0632     if (!tsof.isValid()) {
0633       LogTrace("l1tOmtfEventPrint") << "CandidateSimMuonMatcher::match:" << __LINE__ << " propagation failed"
0634                                     << std::endl;
0635       MatchingResult result;
0636       result.result = MatchingResult::ResultType::propagationFailed;
0637       continue;  //no sense to do matching
0638     }
0639 
0640     /* TODO fix if filling of the deltaPhiPropCandMean and deltaPhiPropCandStdDev is needed
0641     double ptGen = trackingParticle.pt();
0642     if(ptGen >= deltaPhiVertexProp->GetXaxis()->GetXmax())
0643       ptGen = deltaPhiVertexProp->GetXaxis()->GetXmax() - 0.01;
0644 
0645     deltaPhiVertexProp->Fill(ptGen, trackingParticle.momentum().phi() - tsof.globalPosition().phi());
0646 */
0647     unsigned int iCand = 0;
0648     for (auto& muonCand : muonCands) {
0649       //dropping very low quality candidates, as they are fakes usually - but it has no sense, then the results are not conclusive then
0650       /*if(muonCand->hwQual() <= 1)
0651         continue; */
0652 
0653       MatchingResult result;
0654       if (tsof.isValid()) {
0655         result = match(muonCand, ghostBustedProcMuons.at(iCand), trackingParticle, tsof);
0656       }
0657       iCand++;
0658 
0659       if (result.result == MatchingResult::ResultType::matched) {
0660         matchingResults.push_back(result);
0661         matched = true;
0662       }
0663     }
0664 
0665     if (!matched) {  //we are adding result also if it there was no matching to any candidate
0666       MatchingResult result(trackingParticle);
0667       matchingResults.push_back(result);
0668       LogTrace("l1tOmtfEventPrint") << __FUNCTION__ << ":" << __LINE__ << " no matching candidate found" << std::endl;
0669     }
0670   }
0671 
0672   return cleanMatching(matchingResults, muonCands, ghostBustedProcMuons);
0673 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!