Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:27:33 

0001 #include "RecoParticleFlow/PFProducer/interface/PFMuonAlgo.h"
0002 #include "DataFormats/MuonReco/interface/Muon.h"
0003 #include "DataFormats/MuonReco/interface/MuonCocktails.h"
0004 #include "DataFormats/MuonReco/interface/MuonSelectors.h"
0005 #include "DataFormats/ParticleFlowReco/interface/PFBlock.h"
0006 #include "DataFormats/ParticleFlowReco/interface/PFBlockElement.h"
0007 #include "DataFormats/ParticleFlowReco/interface/PFBlockElementGsfTrack.h"
0008 #include "DataFormats/ParticleFlowReco/interface/PFBlockElementTrack.h"
0009 #include "DataFormats/TrackReco/interface/Track.h"
0010 #include <iostream>
0011 #include <memory>
0012 
0013 using namespace std;
0014 using namespace reco;
0015 using namespace boost;
0016 
0017 PFMuonAlgo::PFMuonAlgo(const edm::ParameterSet& iConfig, bool postMuonCleaning)
0018 
0019     : pfCosmicsMuonCleanedCandidates_(std::make_unique<reco::PFCandidateCollection>()),
0020       pfCleanedTrackerAndGlobalMuonCandidates_(std::make_unique<reco::PFCandidateCollection>()),
0021       pfFakeMuonCleanedCandidates_(std::make_unique<reco::PFCandidateCollection>()),
0022       pfPunchThroughMuonCleanedCandidates_(std::make_unique<reco::PFCandidateCollection>()),
0023       pfPunchThroughHadronCleanedCandidates_(std::make_unique<reco::PFCandidateCollection>()),
0024       pfAddedMuonCandidates_(std::make_unique<reco::PFCandidateCollection>()),
0025 
0026       maxDPtOPt_(iConfig.getParameter<double>("maxDPtOPt")),
0027       trackQuality_(reco::TrackBase::qualityByName(iConfig.getParameter<std::string>("trackQuality"))),
0028       errorCompScale_(iConfig.getParameter<double>("ptErrorScale")),
0029 
0030       postCleaning_(postMuonCleaning),  // disable by default (for HLT)
0031       eventFractionCleaning_(iConfig.getParameter<double>("eventFractionForCleaning")),
0032       minPostCleaningPt_(iConfig.getParameter<double>("minPtForPostCleaning")),
0033       eventFactorCosmics_(iConfig.getParameter<double>("eventFactorForCosmics")),
0034       metSigForCleaning_(iConfig.getParameter<double>("metSignificanceForCleaning")),
0035       metSigForRejection_(iConfig.getParameter<double>("metSignificanceForRejection")),
0036       metFactorCleaning_(iConfig.getParameter<double>("metFactorForCleaning")),
0037       eventFractionRejection_(iConfig.getParameter<double>("eventFractionForRejection")),
0038       metFactorRejection_(iConfig.getParameter<double>("metFactorForRejection")),
0039       metFactorHighEta_(iConfig.getParameter<double>("metFactorForHighEta")),
0040       ptFactorHighEta_(iConfig.getParameter<double>("ptFactorForHighEta")),
0041       metFactorFake_(iConfig.getParameter<double>("metFactorForFakes")),
0042       minPunchThroughMomentum_(iConfig.getParameter<double>("minMomentumForPunchThrough")),
0043       minPunchThroughEnergy_(iConfig.getParameter<double>("minEnergyForPunchThrough")),
0044       punchThroughFactor_(iConfig.getParameter<double>("punchThroughFactor")),
0045       punchThroughMETFactor_(iConfig.getParameter<double>("punchThroughMETFactor")),
0046       cosmicRejDistance_(iConfig.getParameter<double>("cosmicRejectionDistance")) {}
0047 
0048 bool PFMuonAlgo::isMuon(const reco::PFBlockElement& elt) {
0049   const auto* eltTrack = dynamic_cast<const reco::PFBlockElementTrack*>(&elt);
0050 
0051   assert(eltTrack);
0052   reco::MuonRef muonRef = eltTrack->muonRef();
0053 
0054   return isMuon(muonRef);
0055 }
0056 
0057 bool PFMuonAlgo::isLooseMuon(const reco::PFBlockElement& elt) {
0058   const reco::PFBlockElementTrack* eltTrack = dynamic_cast<const reco::PFBlockElementTrack*>(&elt);
0059 
0060   assert(eltTrack);
0061 
0062   reco::MuonRef muonRef = eltTrack->muonRef();
0063 
0064   return isLooseMuon(muonRef);
0065 }
0066 
0067 bool PFMuonAlgo::isGlobalTightMuon(const reco::PFBlockElement& elt) {
0068   const reco::PFBlockElementTrack* eltTrack = dynamic_cast<const reco::PFBlockElementTrack*>(&elt);
0069 
0070   assert(eltTrack);
0071   reco::MuonRef muonRef = eltTrack->muonRef();
0072 
0073   return isGlobalTightMuon(muonRef);
0074 }
0075 
0076 bool PFMuonAlgo::isGlobalLooseMuon(const reco::PFBlockElement& elt) {
0077   const reco::PFBlockElementTrack* eltTrack = dynamic_cast<const reco::PFBlockElementTrack*>(&elt);
0078 
0079   assert(eltTrack);
0080   reco::MuonRef muonRef = eltTrack->muonRef();
0081 
0082   return isGlobalLooseMuon(muonRef);
0083 }
0084 
0085 bool PFMuonAlgo::isTrackerTightMuon(const reco::PFBlockElement& elt) {
0086   const reco::PFBlockElementTrack* eltTrack = dynamic_cast<const reco::PFBlockElementTrack*>(&elt);
0087 
0088   assert(eltTrack);
0089   reco::MuonRef muonRef = eltTrack->muonRef();
0090 
0091   return isTrackerTightMuon(muonRef);
0092 }
0093 
0094 bool PFMuonAlgo::isIsolatedMuon(const reco::PFBlockElement& elt) {
0095   const reco::PFBlockElementTrack* eltTrack = dynamic_cast<const reco::PFBlockElementTrack*>(&elt);
0096 
0097   assert(eltTrack);
0098   reco::MuonRef muonRef = eltTrack->muonRef();
0099 
0100   return isIsolatedMuon(muonRef);
0101 }
0102 
0103 bool PFMuonAlgo::isMuon(const reco::MuonRef& muonRef) {
0104   return isGlobalTightMuon(muonRef) || isTrackerTightMuon(muonRef) || isIsolatedMuon(muonRef);
0105 }
0106 
0107 bool PFMuonAlgo::isLooseMuon(const reco::MuonRef& muonRef) {
0108   return (isGlobalLooseMuon(muonRef) || isTrackerLooseMuon(muonRef));
0109 }
0110 
0111 bool PFMuonAlgo::isGlobalTightMuon(const reco::MuonRef& muonRef) {
0112   if (!muonRef.isNonnull())
0113     return false;
0114 
0115   if (!muonRef->isGlobalMuon())
0116     return false;
0117   if (!muonRef->isStandAloneMuon())
0118     return false;
0119 
0120   if (muonRef->isTrackerMuon()) {
0121     bool result = muon::isGoodMuon(*muonRef, muon::GlobalMuonPromptTight);
0122 
0123     bool isTM2DCompatibilityTight = muon::isGoodMuon(*muonRef, muon::TM2DCompatibilityTight);
0124     int nMatches = muonRef->numberOfMatches();
0125     bool quality = nMatches > 2 || isTM2DCompatibilityTight;
0126 
0127     return result && quality;
0128 
0129   } else {
0130     reco::TrackRef standAloneMu = muonRef->standAloneMuon();
0131 
0132     // No tracker muon -> Request a perfect stand-alone muon, or an even better global muon
0133     bool result = false;
0134 
0135     // Check the quality of the stand-alone muon :
0136     // good chi**2 and large number of hits and good pt error
0137     if ((standAloneMu->hitPattern().numberOfValidMuonDTHits() < 22 &&
0138          standAloneMu->hitPattern().numberOfValidMuonCSCHits() < 15) ||
0139         standAloneMu->normalizedChi2() > 10. || standAloneMu->ptError() / standAloneMu->pt() > 0.20) {
0140       result = false;
0141     } else {
0142       reco::TrackRef combinedMu = muonRef->combinedMuon();
0143       reco::TrackRef trackerMu = muonRef->track();
0144 
0145       // If the stand-alone muon is good, check the global muon
0146       if (combinedMu->normalizedChi2() > standAloneMu->normalizedChi2()) {
0147         // If the combined muon is worse than the stand-alone, it
0148         // means that either the corresponding tracker track was not
0149         // reconstructed, or that the sta muon comes from a late
0150         // pion decay (hence with a momentum smaller than the track)
0151         // Take the stand-alone muon only if its momentum is larger
0152         // than that of the track
0153         result = standAloneMu->pt() > trackerMu->pt();
0154       } else {
0155         // If the combined muon is better (and good enough), take the
0156         // global muon
0157         result =
0158             combinedMu->ptError() / combinedMu->pt() < std::min(0.20, standAloneMu->ptError() / standAloneMu->pt());
0159       }
0160     }
0161 
0162     return result;
0163   }
0164 
0165   return false;
0166 }
0167 
0168 bool PFMuonAlgo::isTrackerTightMuon(const reco::MuonRef& muonRef) {
0169   if (!muonRef.isNonnull())
0170     return false;
0171 
0172   if (!muonRef->isTrackerMuon())
0173     return false;
0174 
0175   reco::TrackRef trackerMu = muonRef->track();
0176   const reco::Track& track = *trackerMu;
0177 
0178   unsigned nTrackerHits = track.hitPattern().numberOfValidTrackerHits();
0179 
0180   if (nTrackerHits <= 12)
0181     return false;
0182 
0183   bool isAllArbitrated = muon::isGoodMuon(*muonRef, muon::AllArbitrated);
0184 
0185   bool isTM2DCompatibilityTight = muon::isGoodMuon(*muonRef, muon::TM2DCompatibilityTight);
0186 
0187   if (!isAllArbitrated || !isTM2DCompatibilityTight)
0188     return false;
0189 
0190   if ((trackerMu->ptError() / trackerMu->pt() > 0.10)) {
0191     //std::cout<<" PT ERROR > 10 % "<< trackerMu->pt() <<std::endl;
0192     return false;
0193   }
0194   return true;
0195 }
0196 
0197 bool PFMuonAlgo::isGlobalLooseMuon(const reco::MuonRef& muonRef) {
0198   if (!muonRef.isNonnull())
0199     return false;
0200   if (!muonRef->isGlobalMuon())
0201     return false;
0202   if (!muonRef->isStandAloneMuon())
0203     return false;
0204 
0205   reco::TrackRef standAloneMu = muonRef->standAloneMuon();
0206   reco::TrackRef combinedMu = muonRef->combinedMuon();
0207   reco::TrackRef trackerMu = muonRef->track();
0208 
0209   unsigned nMuonHits =
0210       standAloneMu->hitPattern().numberOfValidMuonDTHits() + 2 * standAloneMu->hitPattern().numberOfValidMuonCSCHits();
0211 
0212   bool quality = false;
0213 
0214   if (muonRef->isTrackerMuon()) {
0215     bool result = combinedMu->normalizedChi2() < 100.;
0216 
0217     bool laststation = muon::isGoodMuon(*muonRef, muon::TMLastStationAngTight);
0218 
0219     int nMatches = muonRef->numberOfMatches();
0220 
0221     quality = laststation && nMuonHits > 12 && nMatches > 1;
0222 
0223     return result && quality;
0224 
0225   } else {
0226     // Check the quality of the stand-alone muon :
0227     // good chi**2 and large number of hits and good pt error
0228     if (nMuonHits <= 15 || standAloneMu->normalizedChi2() > 10. ||
0229         standAloneMu->ptError() / standAloneMu->pt() > 0.20) {
0230       quality = false;
0231     } else {
0232       // If the stand-alone muon is good, check the global muon
0233       if (combinedMu->normalizedChi2() > standAloneMu->normalizedChi2()) {
0234         // If the combined muon is worse than the stand-alone, it
0235         // means that either the corresponding tracker track was not
0236         // reconstructed, or that the sta muon comes from a late
0237         // pion decay (hence with a momentum smaller than the track)
0238         // Take the stand-alone muon only if its momentum is larger
0239         // than that of the track
0240 
0241         // Note that here we even take the standAlone if it has a smaller pT, in contrast to GlobalTight
0242         if (standAloneMu->pt() > trackerMu->pt() || combinedMu->normalizedChi2() < 5.)
0243           quality = true;
0244       } else {
0245         // If the combined muon is better (and good enough), take the
0246         // global muon
0247         if (combinedMu->ptError() / combinedMu->pt() < std::min(0.20, standAloneMu->ptError() / standAloneMu->pt()))
0248           quality = true;
0249       }
0250     }
0251   }
0252 
0253   return quality;
0254 }
0255 
0256 bool PFMuonAlgo::isTrackerLooseMuon(const reco::MuonRef& muonRef) {
0257   if (!muonRef.isNonnull())
0258     return false;
0259   if (!muonRef->isTrackerMuon())
0260     return false;
0261 
0262   reco::TrackRef trackerMu = muonRef->track();
0263 
0264   if (trackerMu->ptError() / trackerMu->pt() > 0.20)
0265     return false;
0266 
0267   // this doesn't seem to be necessary on the small samples looked at, but keep it around as insurance
0268   if (trackerMu->pt() > 20.)
0269     return false;
0270 
0271   bool isAllArbitrated = muon::isGoodMuon(*muonRef, muon::AllArbitrated);
0272   bool isTMLastStationAngTight = muon::isGoodMuon(*muonRef, muon::TMLastStationAngTight);
0273 
0274   bool quality = isAllArbitrated && isTMLastStationAngTight;
0275 
0276   return quality;
0277 }
0278 
0279 bool PFMuonAlgo::isIsolatedMuon(const reco::MuonRef& muonRef) {
0280   if (!muonRef.isNonnull())
0281     return false;
0282   if (!muonRef->isIsolationValid())
0283     return false;
0284 
0285   // Isolated Muons which are missed by standard cuts are nearly always global+tracker
0286   if (!muonRef->isGlobalMuon())
0287     return false;
0288 
0289   // If it's not a tracker muon, only take it if there are valid muon hits
0290 
0291   reco::TrackRef standAloneMu = muonRef->standAloneMuon();
0292 
0293   if (!muonRef->isTrackerMuon()) {
0294     if (standAloneMu->hitPattern().numberOfValidMuonDTHits() == 0 &&
0295         standAloneMu->hitPattern().numberOfValidMuonCSCHits() == 0)
0296       return false;
0297   }
0298 
0299   // for isolation, take the smallest pt available to reject fakes
0300 
0301   reco::TrackRef combinedMu = muonRef->combinedMuon();
0302   double smallestMuPt = combinedMu->pt();
0303 
0304   if (standAloneMu->pt() < smallestMuPt)
0305     smallestMuPt = standAloneMu->pt();
0306 
0307   if (muonRef->isTrackerMuon()) {
0308     reco::TrackRef trackerMu = muonRef->track();
0309     if (trackerMu->pt() < smallestMuPt)
0310       smallestMuPt = trackerMu->pt();
0311   }
0312 
0313   double sumPtR03 = muonRef->isolationR03().sumPt;
0314 
0315   double relIso = sumPtR03 / smallestMuPt;
0316 
0317   if (relIso < 0.1)
0318     return true;
0319   else
0320     return false;
0321 }
0322 
0323 bool PFMuonAlgo::isTightMuonPOG(const reco::MuonRef& muonRef) {
0324   if (!muon::isGoodMuon(*muonRef, muon::GlobalMuonPromptTight))
0325     return false;
0326 
0327   if (!muonRef->isTrackerMuon())
0328     return false;
0329 
0330   if (muonRef->numberOfMatches() < 2)
0331     return false;
0332 
0333   //const reco::TrackRef& combinedMuon = muonRef->combinedMuon();
0334   const reco::TrackRef& combinedMuon = muonRef->globalTrack();
0335 
0336   if (combinedMuon->hitPattern().numberOfValidTrackerHits() < 11)
0337     return false;
0338 
0339   if (combinedMuon->hitPattern().numberOfValidPixelHits() == 0)
0340     return false;
0341 
0342   if (combinedMuon->hitPattern().numberOfValidMuonHits() == 0)
0343     return false;
0344 
0345   return true;
0346 }
0347 
0348 bool PFMuonAlgo::hasValidTrack(const reco::MuonRef& muonRef, bool loose, double maxDPtOPt) {
0349   if (loose)
0350     return !muonTracks(muonRef).empty();
0351   else
0352     return !muonTracks(muonRef, maxDPtOPt).empty();
0353 }
0354 
0355 void PFMuonAlgo::printMuonProperties(const reco::MuonRef& muonRef) {
0356   if (!muonRef.isNonnull())
0357     return;
0358 
0359   bool isGL = muonRef->isGlobalMuon();
0360   bool isTR = muonRef->isTrackerMuon();
0361   bool isST = muonRef->isStandAloneMuon();
0362 
0363   std::cout << " GL: " << isGL << " TR: " << isTR << " ST: " << isST << std::endl;
0364   std::cout << " nMatches " << muonRef->numberOfMatches() << std::endl;
0365 
0366   if (muonRef->isGlobalMuon()) {
0367     reco::TrackRef combinedMu = muonRef->combinedMuon();
0368     std::cout << " GL,  pt: " << combinedMu->pt() << " +/- " << combinedMu->ptError() / combinedMu->pt()
0369               << " chi**2 GBL : " << combinedMu->normalizedChi2() << std::endl;
0370     std::cout << " Total Muon Hits : " << combinedMu->hitPattern().numberOfValidMuonHits() << "/"
0371               << combinedMu->hitPattern().numberOfLostMuonHits()
0372               << " DT Hits : " << combinedMu->hitPattern().numberOfValidMuonDTHits() << "/"
0373               << combinedMu->hitPattern().numberOfLostMuonDTHits()
0374               << " CSC Hits : " << combinedMu->hitPattern().numberOfValidMuonCSCHits() << "/"
0375               << combinedMu->hitPattern().numberOfLostMuonCSCHits()
0376               << " RPC Hits : " << combinedMu->hitPattern().numberOfValidMuonRPCHits() << "/"
0377               << combinedMu->hitPattern().numberOfLostMuonRPCHits() << std::endl;
0378 
0379     std::cout << "  # of Valid Tracker Hits " << combinedMu->hitPattern().numberOfValidTrackerHits() << std::endl;
0380     std::cout << "  # of Valid Pixel Hits " << combinedMu->hitPattern().numberOfValidPixelHits() << std::endl;
0381   }
0382   if (muonRef->isStandAloneMuon()) {
0383     reco::TrackRef standAloneMu = muonRef->standAloneMuon();
0384     std::cout << " ST,  pt: " << standAloneMu->pt() << " +/- " << standAloneMu->ptError() / standAloneMu->pt()
0385               << " eta : " << standAloneMu->eta()
0386               << " DT Hits : " << standAloneMu->hitPattern().numberOfValidMuonDTHits() << "/"
0387               << standAloneMu->hitPattern().numberOfLostMuonDTHits()
0388               << " CSC Hits : " << standAloneMu->hitPattern().numberOfValidMuonCSCHits() << "/"
0389               << standAloneMu->hitPattern().numberOfLostMuonCSCHits()
0390               << " RPC Hits : " << standAloneMu->hitPattern().numberOfValidMuonRPCHits() << "/"
0391               << standAloneMu->hitPattern().numberOfLostMuonRPCHits()
0392               << " chi**2 STA : " << standAloneMu->normalizedChi2() << std::endl;
0393   }
0394 
0395   if (muonRef->isTrackerMuon()) {
0396     reco::TrackRef trackerMu = muonRef->track();
0397     const reco::Track& track = *trackerMu;
0398     std::cout << " TR,  pt: " << trackerMu->pt() << " +/- " << trackerMu->ptError() / trackerMu->pt()
0399               << " chi**2 TR : " << trackerMu->normalizedChi2() << std::endl;
0400     std::cout << " nTrackerHits " << track.hitPattern().numberOfValidTrackerHits() << std::endl;
0401     std::cout << "TMLastStationAngLoose               " << muon::isGoodMuon(*muonRef, muon::TMLastStationAngLoose)
0402               << std::endl
0403               << "TMLastStationAngTight               " << muon::isGoodMuon(*muonRef, muon::TMLastStationAngTight)
0404               << std::endl
0405               << "TMLastStationLoose               " << muon::isGoodMuon(*muonRef, muon::TMLastStationLoose)
0406               << std::endl
0407               << "TMLastStationTight               " << muon::isGoodMuon(*muonRef, muon::TMLastStationTight)
0408               << std::endl
0409               << "TMOneStationLoose                " << muon::isGoodMuon(*muonRef, muon::TMOneStationLoose) << std::endl
0410               << "TMOneStationTight                " << muon::isGoodMuon(*muonRef, muon::TMOneStationTight) << std::endl
0411               << "TMLastStationOptimizedLowPtLoose "
0412               << muon::isGoodMuon(*muonRef, muon::TMLastStationOptimizedLowPtLoose) << std::endl
0413               << "TMLastStationOptimizedLowPtTight "
0414               << muon::isGoodMuon(*muonRef, muon::TMLastStationOptimizedLowPtTight) << std::endl
0415               << "TMLastStationOptimizedBarrelLowPtLoose "
0416               << muon::isGoodMuon(*muonRef, muon::TMLastStationOptimizedBarrelLowPtLoose) << std::endl
0417               << "TMLastStationOptimizedBarrelLowPtTight "
0418               << muon::isGoodMuon(*muonRef, muon::TMLastStationOptimizedBarrelLowPtTight) << std::endl
0419               << std::endl;
0420   }
0421 
0422   std::cout << "TM2DCompatibilityLoose           " << muon::isGoodMuon(*muonRef, muon::TM2DCompatibilityLoose)
0423             << std::endl
0424             << "TM2DCompatibilityTight           " << muon::isGoodMuon(*muonRef, muon::TM2DCompatibilityTight)
0425             << std::endl;
0426 
0427   if (muonRef->isGlobalMuon() && muonRef->isTrackerMuon() && muonRef->isStandAloneMuon()) {
0428     reco::TrackRef combinedMu = muonRef->combinedMuon();
0429     reco::TrackRef trackerMu = muonRef->track();
0430     reco::TrackRef standAloneMu = muonRef->standAloneMuon();
0431 
0432     double sigmaCombined = combinedMu->ptError() / (combinedMu->pt() * combinedMu->pt());
0433     double sigmaTracker = trackerMu->ptError() / (trackerMu->pt() * trackerMu->pt());
0434     double sigmaStandAlone = standAloneMu->ptError() / (standAloneMu->pt() * standAloneMu->pt());
0435 
0436     bool combined = combinedMu->ptError() / combinedMu->pt() < 0.20;
0437     bool tracker = trackerMu->ptError() / trackerMu->pt() < 0.20;
0438     bool standAlone = standAloneMu->ptError() / standAloneMu->pt() < 0.20;
0439 
0440     double delta1 = combined && tracker ? fabs(1. / combinedMu->pt() - 1. / trackerMu->pt()) /
0441                                               sqrt(sigmaCombined * sigmaCombined + sigmaTracker * sigmaTracker)
0442                                         : 100.;
0443     double delta2 = combined && standAlone ? fabs(1. / combinedMu->pt() - 1. / standAloneMu->pt()) /
0444                                                  sqrt(sigmaCombined * sigmaCombined + sigmaStandAlone * sigmaStandAlone)
0445                                            : 100.;
0446     double delta3 = standAlone && tracker ? fabs(1. / standAloneMu->pt() - 1. / trackerMu->pt()) /
0447                                                 sqrt(sigmaStandAlone * sigmaStandAlone + sigmaTracker * sigmaTracker)
0448                                           : 100.;
0449 
0450     double delta =
0451         standAloneMu->hitPattern().numberOfValidMuonDTHits() + standAloneMu->hitPattern().numberOfValidMuonCSCHits() > 0
0452             ? std::min(delta3, std::min(delta1, delta2))
0453             : std::max(delta3, std::max(delta1, delta2));
0454 
0455     std::cout << "delta = " << delta << " delta1 " << delta1 << " delta2 " << delta2 << " delta3 " << delta3
0456               << std::endl;
0457 
0458     double ratio = combinedMu->ptError() / combinedMu->pt() / (trackerMu->ptError() / trackerMu->pt());
0459     //if ( ratio > 2. && delta < 3. ) std::cout << "ALARM ! " << ratio << ", " << delta << std::endl;
0460     std::cout << " ratio " << ratio << " combined mu pt " << combinedMu->pt() << std::endl;
0461     //bool quality3 =  ( combinedMu->pt() < 50. || ratio < 2. ) && delta <  3.;
0462   }
0463 
0464   double sumPtR03 = muonRef->isolationR03().sumPt;
0465   double emEtR03 = muonRef->isolationR03().emEt;
0466   double hadEtR03 = muonRef->isolationR03().hadEt;
0467   double relIsoR03 = (sumPtR03 + emEtR03 + hadEtR03) / muonRef->pt();
0468   double sumPtR05 = muonRef->isolationR05().sumPt;
0469   double emEtR05 = muonRef->isolationR05().emEt;
0470   double hadEtR05 = muonRef->isolationR05().hadEt;
0471   double relIsoR05 = (sumPtR05 + emEtR05 + hadEtR05) / muonRef->pt();
0472   std::cout << " 0.3 Radion Rel Iso: " << relIsoR03 << " sumPt " << sumPtR03 << " emEt " << emEtR03 << " hadEt "
0473             << hadEtR03 << std::endl;
0474   std::cout << " 0.5 Radion Rel Iso: " << relIsoR05 << " sumPt " << sumPtR05 << " emEt " << emEtR05 << " hadEt "
0475             << hadEtR05 << std::endl;
0476   return;
0477 }
0478 
0479 int PFMuonAlgo::muAssocToTrack(const reco::TrackRef& trackref, const reco::MuonCollection& muons) {
0480   auto muon = std::find_if(muons.cbegin(), muons.cend(), [&](const reco::Muon& m) {
0481     return (m.track().isNonnull() && m.track() == trackref);
0482   });
0483   return (muon != muons.cend() ? std::distance(muons.cbegin(), muon) : -1);
0484 }
0485 
0486 std::vector<reco::Muon::MuonTrackTypePair> PFMuonAlgo::muonTracks(const reco::MuonRef& muon,
0487                                                                   double maxDPtOPt,
0488                                                                   bool includeSA) {
0489   std::vector<reco::Muon::MuonTrackTypePair> out;
0490 
0491   if (muon->globalTrack().isNonnull() && muon->globalTrack()->pt() > 0)
0492     if (muon->globalTrack()->ptError() / muon->globalTrack()->pt() < maxDPtOPt)
0493       out.emplace_back(muon->globalTrack(), reco::Muon::CombinedTrack);
0494 
0495   if (muon->innerTrack().isNonnull() && muon->innerTrack()->pt() > 0)
0496     if (muon->innerTrack()->ptError() / muon->innerTrack()->pt() < maxDPtOPt)  //Here Loose!@
0497       out.emplace_back(muon->innerTrack(), reco::Muon::InnerTrack);
0498 
0499   bool pickyExists = false;
0500   double pickyDpt = 99999.;
0501   if (muon->pickyTrack().isNonnull() && muon->pickyTrack()->pt() > 0) {
0502     pickyDpt = muon->pickyTrack()->ptError() / muon->pickyTrack()->pt();
0503     if (pickyDpt < maxDPtOPt)
0504       out.emplace_back(muon->pickyTrack(), reco::Muon::Picky);
0505     pickyExists = true;
0506   }
0507 
0508   bool dytExists = false;
0509   double dytDpt = 99999.;
0510   if (muon->dytTrack().isNonnull() && muon->dytTrack()->pt() > 0) {
0511     dytDpt = muon->dytTrack()->ptError() / muon->dytTrack()->pt();
0512     if (dytDpt < maxDPtOPt)
0513       out.emplace_back(muon->dytTrack(), reco::Muon::DYT);
0514     dytExists = true;
0515   }
0516 
0517   //Magic: TPFMS is not a really good track especially under misalignment
0518   //IT is kind of crap because if mu system is displaced it can make a change
0519   //So allow TPFMS if there is no picky or the error of tpfms is better than picky
0520   //AND if there is no DYT or the error of tpfms is better than DYT
0521   if (muon->tpfmsTrack().isNonnull() && muon->tpfmsTrack()->pt() > 0) {
0522     double tpfmsDpt = muon->tpfmsTrack()->ptError() / muon->tpfmsTrack()->pt();
0523     if (((pickyExists && tpfmsDpt < pickyDpt) || (!pickyExists)) &&
0524         ((dytExists && tpfmsDpt < dytDpt) || (!dytExists)) && tpfmsDpt < maxDPtOPt)
0525       out.emplace_back(muon->tpfmsTrack(), reco::Muon::TPFMS);
0526   }
0527 
0528   if (includeSA && muon->outerTrack().isNonnull())
0529     if (muon->outerTrack()->ptError() / muon->outerTrack()->pt() < maxDPtOPt)
0530       out.emplace_back(muon->outerTrack(), reco::Muon::OuterTrack);
0531 
0532   return out;
0533 }
0534 
0535 ////////////////////////////////////////////////////////////////////////////////////
0536 
0537 bool PFMuonAlgo::reconstructMuon(reco::PFCandidate& candidate, const reco::MuonRef& muon, bool allowLoose) {
0538   using namespace std;
0539   using namespace reco;
0540 
0541   if (!muon.isNonnull())
0542     return false;
0543 
0544   bool isMu = false;
0545 
0546   if (allowLoose)
0547     isMu = isMuon(muon) || isLooseMuon(muon);
0548   else
0549     isMu = isMuon(muon);
0550 
0551   if (!isMu)
0552     return false;
0553 
0554   //get the valid tracks(without standalone except we allow loose muons)
0555   //MIKE: Here we need to be careful. If we have a muon inside a dense
0556   //jet environment often the track is badly measured. In this case
0557   //we should not apply Dpt/Pt<1
0558 
0559   std::vector<reco::Muon::MuonTrackTypePair> validTracks = muonTracks(muon, maxDPtOPt_);
0560   if (!allowLoose)
0561     validTracks = muonTracks(muon, maxDPtOPt_);
0562   else
0563     validTracks = muonTracks(muon);
0564 
0565   if (validTracks.empty())
0566     return false;
0567 
0568   //check what is the track used.Rerun TuneP
0569   reco::Muon::MuonTrackTypePair bestTrackPair = muon::tevOptimized(*muon);
0570 
0571   TrackRef bestTrack = bestTrackPair.first;
0572   MuonTrackType trackType = bestTrackPair.second;
0573 
0574   MuonTrackTypePair trackPairWithSmallestError = getTrackWithSmallestError(validTracks);
0575   TrackRef trackWithSmallestError = trackPairWithSmallestError.first;
0576 
0577   if (trackType == reco::Muon::InnerTrack &&
0578       (!bestTrack->quality(trackQuality_) ||
0579        bestTrack->ptError() / bestTrack->pt() >
0580            errorCompScale_ * trackWithSmallestError->ptError() / trackWithSmallestError->pt())) {
0581     bestTrack = trackWithSmallestError;
0582     trackType = trackPairWithSmallestError.second;
0583   } else if (trackType != reco::Muon::InnerTrack &&
0584              bestTrack->ptError() / bestTrack->pt() >
0585                  errorCompScale_ * trackWithSmallestError->ptError() / trackWithSmallestError->pt()) {
0586     bestTrack = trackWithSmallestError;
0587     trackType = trackPairWithSmallestError.second;
0588   }
0589 
0590   changeTrack(candidate, std::make_pair(bestTrack, trackType));
0591   candidate.setMuonRef(muon);
0592 
0593   return true;
0594 }
0595 
0596 void PFMuonAlgo::changeTrack(reco::PFCandidate& candidate, const MuonTrackTypePair& track) {
0597   using namespace reco;
0598   reco::TrackRef bestTrack = track.first;
0599   MuonTrackType trackType = track.second;
0600   //OK Now redefine the canddiate with that track
0601   double px = bestTrack->px();
0602   double py = bestTrack->py();
0603   double pz = bestTrack->pz();
0604   double energy = sqrt(bestTrack->p() * bestTrack->p() + 0.1057 * 0.1057);
0605 
0606   candidate.setCharge(bestTrack->charge() > 0 ? 1 : -1);
0607   candidate.setP4(math::XYZTLorentzVector(px, py, pz, energy));
0608   candidate.setParticleType(reco::PFCandidate::mu);
0609   //    candidate.setTrackRef( bestTrack );
0610   candidate.setMuonTrackType(trackType);
0611   candidate.setVertex(bestTrack->vertex());
0612 }
0613 
0614 reco::Muon::MuonTrackTypePair PFMuonAlgo::getTrackWithSmallestError(
0615     const std::vector<reco::Muon::MuonTrackTypePair>& tracks) {
0616   TrackPtErrorSorter sorter;
0617   return *std::min_element(tracks.begin(), tracks.end(), sorter);
0618 }
0619 
0620 void PFMuonAlgo::estimateEventQuantities(const reco::PFCandidateCollection* pfc) {
0621   //SUM ET
0622   sumetPU_ = 0.0;
0623   METX_ = 0.;
0624   METY_ = 0.;
0625   sumet_ = 0.0;
0626   for (reco::PFCandidateCollection::const_iterator i = pfc->begin(); i != pfc->end(); ++i) {
0627     sumet_ += i->pt();
0628     METX_ += i->px();
0629     METY_ += i->py();
0630   }
0631 }
0632 
0633 void PFMuonAlgo::postClean(reco::PFCandidateCollection* cands) {
0634   using namespace std;
0635   using namespace reco;
0636   if (!postCleaning_)
0637     return;
0638 
0639   //Initialize vectors
0640 
0641   if (pfCosmicsMuonCleanedCandidates_.get())
0642     pfCosmicsMuonCleanedCandidates_->clear();
0643   else
0644     pfCosmicsMuonCleanedCandidates_ = std::make_unique<reco::PFCandidateCollection>();
0645 
0646   if (pfCleanedTrackerAndGlobalMuonCandidates_.get())
0647     pfCleanedTrackerAndGlobalMuonCandidates_->clear();
0648   else
0649     pfCleanedTrackerAndGlobalMuonCandidates_ = std::make_unique<reco::PFCandidateCollection>();
0650 
0651   if (pfFakeMuonCleanedCandidates_.get())
0652     pfFakeMuonCleanedCandidates_->clear();
0653   else
0654     pfFakeMuonCleanedCandidates_ = std::make_unique<reco::PFCandidateCollection>();
0655 
0656   if (pfPunchThroughMuonCleanedCandidates_.get())
0657     pfPunchThroughMuonCleanedCandidates_->clear();
0658   else
0659     pfPunchThroughMuonCleanedCandidates_ = std::make_unique<reco::PFCandidateCollection>();
0660 
0661   if (pfPunchThroughHadronCleanedCandidates_.get())
0662     pfPunchThroughHadronCleanedCandidates_->clear();
0663   else
0664     pfPunchThroughHadronCleanedCandidates_ = std::make_unique<reco::PFCandidateCollection>();
0665 
0666   pfPunchThroughHadronCleanedCandidates_->clear();
0667 
0668   maskedIndices_.clear();
0669 
0670   //Estimate MET and SumET
0671 
0672   estimateEventQuantities(cands);
0673 
0674   std::vector<int> muons;
0675   std::vector<int> cosmics;
0676   //get the muons
0677   for (unsigned int i = 0; i < cands->size(); ++i) {
0678     const reco::PFCandidate& cand = (*cands)[i];
0679     if (cand.particleId() == reco::PFCandidate::mu)
0680       muons.push_back(i);
0681   }
0682   //Then sort the muon indicess by decsending pt
0683 
0684   IndexPtComparator comparator(cands);
0685   std::sort(muons.begin(), muons.end(), comparator);
0686 
0687   //first kill cosmics
0688   double METXCosmics = 0;
0689   double METYCosmics = 0;
0690   double SUMETCosmics = 0.0;
0691 
0692   for (unsigned int i = 0; i < muons.size(); ++i) {
0693     const PFCandidate& pfc = cands->at(muons[i]);
0694     double origin = 0.0;
0695     if (!vertices_->empty() && vertices_->at(0).isValid() && !vertices_->at(0).isFake())
0696       origin = pfc.muonRef()->muonBestTrack()->dxy(vertices_->at(0).position());
0697 
0698     if (origin > cosmicRejDistance_) {
0699       cosmics.push_back(muons[i]);
0700       METXCosmics += pfc.px();
0701       METYCosmics += pfc.py();
0702       SUMETCosmics += pfc.pt();
0703     }
0704   }
0705   double MET2Cosmics = METXCosmics * METXCosmics + METYCosmics * METYCosmics;
0706 
0707   if (SUMETCosmics > (sumet_ - sumetPU_) / eventFactorCosmics_ && MET2Cosmics < METX_ * METX_ + METY_ * METY_)
0708     for (unsigned int i = 0; i < cosmics.size(); ++i) {
0709       maskedIndices_.push_back(cosmics[i]);
0710       pfCosmicsMuonCleanedCandidates_->push_back(cands->at(cosmics[i]));
0711     }
0712 
0713   //Loop on the muons candidates and clean
0714   for (unsigned int i = 0; i < muons.size(); ++i) {
0715     if (cleanMismeasured(cands->at(muons[i]), muons[i]))
0716       continue;
0717     cleanPunchThroughAndFakes(cands->at(muons[i]), cands, muons[i]);
0718   }
0719 
0720   //OK Now do the hard job ->remove the candidates that were cleaned
0721   removeDeadCandidates(cands, maskedIndices_);
0722 }
0723 
0724 void PFMuonAlgo::addMissingMuons(edm::Handle<reco::MuonCollection> muons, reco::PFCandidateCollection* cands) {
0725   if (!postCleaning_)
0726     return;
0727 
0728   if (pfAddedMuonCandidates_.get())
0729     pfAddedMuonCandidates_->clear();
0730   else
0731     pfAddedMuonCandidates_ = std::make_unique<reco::PFCandidateCollection>();
0732 
0733   // do not recover muons if they are in the veto list
0734   typedef std::pair<edm::ProductID, unsigned> TrackProdIDKey;
0735   std::vector<TrackProdIDKey> vetoed;
0736   if (vetoes_) {
0737     for (const auto& veto : *vetoes_) {
0738       if (veto.trackRef().isNull())
0739         continue;
0740       vetoed.emplace_back(veto.trackRef().id(), veto.trackRef().key());
0741     }
0742     std::sort(vetoed.begin(), vetoed.end());
0743   }
0744 
0745   for (unsigned imu = 0; imu < muons->size(); ++imu) {
0746     reco::MuonRef muonRef(muons, imu);
0747     bool used = false;
0748     bool hadron = false;
0749     for (unsigned i = 0; i < cands->size(); i++) {
0750       const PFCandidate& pfc = cands->at(i);
0751       if (!pfc.trackRef().isNonnull())
0752         continue;
0753       if (pfc.trackRef().isNonnull() && pfc.trackRef() == muonRef->track())
0754         hadron = true;
0755       if (!pfc.muonRef().isNonnull())
0756         continue;
0757 
0758       if (pfc.muonRef()->innerTrack() == muonRef->innerTrack())
0759         used = true;
0760       else {
0761         // Check if the stand-alone muon is not a spurious copy of an existing muon
0762         // (Protection needed for HLT)
0763         if (pfc.muonRef()->isStandAloneMuon() && muonRef->isStandAloneMuon()) {
0764           double dEta = pfc.muonRef()->standAloneMuon()->eta() - muonRef->standAloneMuon()->eta();
0765           double dPhi = pfc.muonRef()->standAloneMuon()->phi() - muonRef->standAloneMuon()->phi();
0766           double dR = sqrt(dEta * dEta + dPhi * dPhi);
0767           if (dR < 0.005) {
0768             used = true;
0769           }
0770         }
0771       }
0772 
0773       if (used)
0774         break;
0775     }
0776 
0777     if (!vetoed.empty()) {
0778       TrackProdIDKey trk = std::make_pair(muonRef->track().id(), muonRef->track().key());
0779       auto lower = std::lower_bound(vetoed.begin(), vetoed.end(), trk);
0780       bool inVetoList = (lower != vetoed.end() && *lower == trk);
0781       if (inVetoList)
0782         used = true;
0783     }
0784 
0785     if (used || hadron || (!muonRef.isNonnull()))
0786       continue;
0787 
0788     TrackMETComparator comparator(METX_, METY_);
0789     //Low pt dont need to be cleaned
0790 
0791     std::vector<reco::Muon::MuonTrackTypePair> tracks = muonTracks(muonRef, maxDPtOPt_, true);
0792     //If there is at least 1 track choice  try to change the track
0793     if (!tracks.empty()) {
0794       //Find tracks that change dramatically MET or Pt
0795       std::vector<reco::Muon::MuonTrackTypePair> tracksThatChangeMET = tracksPointingAtMET(tracks);
0796       //From those tracks get the one with smallest MET
0797       if (!tracksThatChangeMET.empty()) {
0798         reco::Muon::MuonTrackTypePair bestTrackType =
0799             *std::min_element(tracksThatChangeMET.begin(), tracksThatChangeMET.end(), comparator);
0800 
0801         //Make sure it is not cosmic
0802         if ((vertices_->empty()) || bestTrackType.first->dz(vertices_->at(0).position()) < cosmicRejDistance_) {
0803           //make a pfcandidate
0804           int charge = bestTrackType.first->charge() > 0 ? 1 : -1;
0805           math::XYZTLorentzVector momentum(bestTrackType.first->px(),
0806                                            bestTrackType.first->py(),
0807                                            bestTrackType.first->pz(),
0808                                            sqrt(bestTrackType.first->p() * bestTrackType.first->p() + 0.1057 * 0.1057));
0809 
0810           cands->push_back(PFCandidate(charge, momentum, reco::PFCandidate::mu));
0811 
0812           changeTrack(cands->back(), bestTrackType);
0813 
0814           if (muonRef->track().isNonnull())
0815             cands->back().setTrackRef(muonRef->track());
0816 
0817           cands->back().setMuonRef(muonRef);
0818 
0819           pfAddedMuonCandidates_->push_back(cands->back());
0820         }
0821       }
0822     }
0823   }
0824 }
0825 
0826 std::pair<double, double> PFMuonAlgo::getMinMaxMET2(const reco::PFCandidate& pfc) {
0827   std::vector<reco::Muon::MuonTrackTypePair> tracks = muonTracks((pfc.muonRef()), maxDPtOPt_, true);
0828 
0829   double METXNO = METX_ - pfc.px();
0830   double METYNO = METY_ - pfc.py();
0831   std::vector<double> met2;
0832   for (unsigned int i = 0; i < tracks.size(); ++i) {
0833     met2.push_back(pow(METXNO + tracks.at(i).first->px(), 2) + pow(METYNO + tracks.at(i).first->py(), 2));
0834   }
0835 
0836   //PROTECT for cases of only one track. If there is only one track it will crash .
0837   //Has never happened but could likely happen!
0838 
0839   if (tracks.size() > 1)
0840     return std::make_pair(*std::min_element(met2.begin(), met2.end()), *std::max_element(met2.begin(), met2.end()));
0841   else
0842     return std::make_pair(0, 0);
0843 }
0844 
0845 bool PFMuonAlgo::cleanMismeasured(reco::PFCandidate& pfc, unsigned int i) {
0846   using namespace std;
0847   using namespace reco;
0848   bool cleaned = false;
0849 
0850   //First define the MET without this guy
0851   double METNOX = METX_ - pfc.px();
0852   double METNOY = METY_ - pfc.py();
0853   double SUMETNO = sumet_ - pfc.pt();
0854 
0855   TrackMETComparator comparator(METNOX, METNOY);
0856   //Low pt dont need to be cleaned
0857   if (pfc.pt() < minPostCleaningPt_)
0858     return false;
0859   std::vector<reco::Muon::MuonTrackTypePair> tracks = muonTracks(pfc.muonRef(), maxDPtOPt_, false);
0860 
0861   //If there is more than 1 track choice  try to change the track
0862   if (tracks.size() > 1) {
0863     //Find tracks that change dramatically MET or Pt
0864     std::vector<reco::Muon::MuonTrackTypePair> tracksThatChangeMET = tracksWithBetterMET(tracks, pfc);
0865     //From those tracks get the one with smallest MET
0866     if (!tracksThatChangeMET.empty()) {
0867       reco::Muon::MuonTrackTypePair bestTrackType =
0868           *std::min_element(tracksThatChangeMET.begin(), tracksThatChangeMET.end(), comparator);
0869       changeTrack(pfc, bestTrackType);
0870 
0871       pfCleanedTrackerAndGlobalMuonCandidates_->push_back(pfc);
0872       //update eventquantities
0873       METX_ = METNOX + pfc.px();
0874       METY_ = METNOY + pfc.py();
0875       sumet_ = SUMETNO + pfc.pt();
0876     }
0877   }
0878 
0879   //Now attempt to kill it
0880   if (!(pfc.muonRef()->isGlobalMuon() && pfc.muonRef()->isTrackerMuon())) {
0881     //define MET significance and SUM ET
0882     double MET2 = METX_ * METX_ + METY_ * METY_;
0883     double newMET2 = METNOX * METNOX + METNOY * METNOY;
0884     double METSig = sqrt(MET2) / sqrt(sumet_ - sumetPU_);
0885     if (METSig > metSigForRejection_)
0886       if ((newMET2 < MET2 / metFactorRejection_) &&
0887           ((SUMETNO - sumetPU_) / (sumet_ - sumetPU_) < eventFractionRejection_)) {
0888         pfFakeMuonCleanedCandidates_->push_back(pfc);
0889         maskedIndices_.push_back(i);
0890         METX_ = METNOX;
0891         METY_ = METNOY;
0892         sumet_ = SUMETNO;
0893         cleaned = true;
0894       }
0895   }
0896   return cleaned;
0897 }
0898 
0899 std::vector<reco::Muon::MuonTrackTypePair> PFMuonAlgo::tracksWithBetterMET(
0900     const std::vector<reco::Muon::MuonTrackTypePair>& tracks, const reco::PFCandidate& pfc) {
0901   std::vector<reco::Muon::MuonTrackTypePair> outputTracks;
0902 
0903   double METNOX = METX_ - pfc.px();
0904   double METNOY = METY_ - pfc.py();
0905   double SUMETNO = sumet_ - pfc.pt();
0906   double MET2 = METX_ * METX_ + METY_ * METY_;
0907   double newMET2 = 0.0;
0908   double newSUMET = 0.0;
0909   double METSIG = sqrt(MET2) / sqrt(sumet_ - sumetPU_);
0910 
0911   if (METSIG > metSigForCleaning_)
0912     for (unsigned int i = 0; i < tracks.size(); ++i) {
0913       //calculate new SUM ET and MET2
0914       newSUMET = SUMETNO + tracks.at(i).first->pt() - sumetPU_;
0915       newMET2 = pow(METNOX + tracks.at(i).first->px(), 2) + pow(METNOY + tracks.at(i).first->py(), 2);
0916 
0917       if (newSUMET / (sumet_ - sumetPU_) > eventFractionCleaning_ && newMET2 < MET2 / metFactorCleaning_)
0918         outputTracks.push_back(tracks.at(i));
0919     }
0920 
0921   return outputTracks;
0922 }
0923 
0924 std::vector<reco::Muon::MuonTrackTypePair> PFMuonAlgo::tracksPointingAtMET(
0925     const std::vector<reco::Muon::MuonTrackTypePair>& tracks) {
0926   std::vector<reco::Muon::MuonTrackTypePair> outputTracks;
0927 
0928   double newMET2 = 0.0;
0929 
0930   for (unsigned int i = 0; i < tracks.size(); ++i) {
0931     //calculate new SUM ET and MET2
0932     newMET2 = pow(METX_ + tracks.at(i).first->px(), 2) + pow(METY_ + tracks.at(i).first->py(), 2);
0933 
0934     if (newMET2 < (METX_ * METX_ + METY_ * METY_) / metFactorCleaning_)
0935       outputTracks.push_back(tracks.at(i));
0936   }
0937 
0938   return outputTracks;
0939 }
0940 
0941 void PFMuonAlgo::setInputsForCleaning(reco::VertexCollection const& vertices) { vertices_ = &vertices; }
0942 
0943 bool PFMuonAlgo::cleanPunchThroughAndFakes(reco::PFCandidate& pfc,
0944                                            reco::PFCandidateCollection* cands,
0945                                            unsigned int imu) {
0946   using namespace reco;
0947 
0948   bool cleaned = false;
0949 
0950   if (pfc.pt() < minPostCleaningPt_)
0951     return false;
0952 
0953   double METXNO = METX_ - pfc.pt();
0954   double METYNO = METY_ - pfc.pt();
0955   double MET2NO = METXNO * METXNO + METYNO * METYNO;
0956   double MET2 = METX_ * METX_ + METY_ * METY_;
0957   bool fake1 = false;
0958 
0959   std::pair<double, double> met2 = getMinMaxMET2(pfc);
0960 
0961   //Check for Fakes at high pseudorapidity
0962   if (pfc.muonRef()->standAloneMuon().isNonnull())
0963     fake1 = fabs(pfc.eta()) > 2.15 && met2.first < met2.second / 2 && MET2NO < MET2 / metFactorHighEta_ &&
0964             pfc.muonRef()->standAloneMuon()->pt() < pfc.pt() / ptFactorHighEta_;
0965 
0966   double factor = std::max(2., 2000. / (sumet_ - pfc.pt() - sumetPU_));
0967   bool fake2 =
0968       (pfc.pt() / (sumet_ - sumetPU_) < 0.25 && MET2NO < MET2 / metFactorFake_ && met2.first < met2.second / factor);
0969 
0970   bool punchthrough = pfc.p() > minPunchThroughMomentum_ && pfc.rawHcalEnergy() > minPunchThroughEnergy_ &&
0971                       pfc.rawEcalEnergy() + pfc.rawHcalEnergy() > pfc.p() / punchThroughFactor_ &&
0972                       !isIsolatedMuon(pfc.muonRef()) && MET2NO < MET2 / punchThroughMETFactor_;
0973 
0974   if (fake1 || fake2 || punchthrough) {
0975     // Find the block of the muon
0976     const PFCandidate::ElementsInBlocks& eleInBlocks = pfc.elementsInBlocks();
0977     if (!eleInBlocks.empty()) {
0978       PFBlockRef blockRefMuon = eleInBlocks[0].first;
0979       unsigned indexMuon = eleInBlocks[0].second;
0980       if (eleInBlocks.size() > 1)
0981         indexMuon = eleInBlocks[1].second;
0982 
0983       // Check if the muon gave rise to a neutral hadron
0984       double iHad = 1E9;
0985       bool hadron = false;
0986       for (unsigned i = imu + 1; i < cands->size(); ++i) {
0987         const PFCandidate& pfcn = cands->at(i);
0988         const PFCandidate::ElementsInBlocks& ele = pfcn.elementsInBlocks();
0989         if (ele.empty()) {
0990           continue;
0991         }
0992         PFBlockRef blockRefHadron = ele[0].first;
0993         unsigned indexHadron = ele[0].second;
0994         // We are out of the block -> exit the loop
0995         if (blockRefHadron.key() != blockRefMuon.key())
0996           break;
0997         // Check that this particle is a neutral hadron
0998         if (indexHadron == indexMuon && pfcn.particleId() == reco::PFCandidate::h0) {
0999           iHad = i;
1000           hadron = true;
1001         }
1002         if (hadron)
1003           break;
1004       }
1005 
1006       if (hadron) {
1007         double rescaleFactor = cands->at(iHad).p() / cands->at(imu).p();
1008         METX_ -= cands->at(imu).px() + cands->at(iHad).px();
1009         METY_ -= cands->at(imu).py() + cands->at(iHad).py();
1010         sumet_ -= cands->at(imu).pt();
1011         cands->at(imu).rescaleMomentum(rescaleFactor);
1012         maskedIndices_.push_back(iHad);
1013         pfPunchThroughHadronCleanedCandidates_->push_back(cands->at(iHad));
1014         cands->at(imu).setParticleType(reco::PFCandidate::h);
1015         pfPunchThroughMuonCleanedCandidates_->push_back(cands->at(imu));
1016         METX_ += cands->at(imu).px();
1017         METY_ += cands->at(imu).py();
1018         sumet_ += cands->at(imu).pt();
1019 
1020       } else if (fake1 || fake2) {
1021         METX_ -= cands->at(imu).px();
1022         METY_ -= cands->at(imu).py();
1023         sumet_ -= cands->at(imu).pt();
1024         maskedIndices_.push_back(imu);
1025         pfFakeMuonCleanedCandidates_->push_back(cands->at(imu));
1026         cleaned = true;
1027       }
1028     }
1029   }
1030   return cleaned;
1031 }
1032 
1033 void PFMuonAlgo::removeDeadCandidates(reco::PFCandidateCollection* obj, const std::vector<unsigned int>& indices) {
1034   size_t N = indices.size();
1035   size_t collSize = obj->size();
1036 
1037   for (size_t i = 0; i < N; ++i)
1038     obj->at(indices.at(i)) = obj->at(collSize - i - 1);
1039 
1040   obj->resize(collSize - indices.size());
1041 }
1042 
1043 void PFMuonAlgo::fillPSetDescription(edm::ParameterSetDescription& iDesc) {
1044   // Muon ID and post cleaning parameters
1045   iDesc.add<double>("maxDPtOPt", 1.0);
1046   iDesc.add<std::string>("trackQuality", "highPurity");
1047   iDesc.add<double>("ptErrorScale", 8.0);
1048 
1049   iDesc.add<double>("eventFractionForCleaning", 0.5);
1050   iDesc.add<double>("minPtForPostCleaning", 20.0);
1051   iDesc.add<double>("eventFactorForCosmics", 10.0);
1052   iDesc.add<double>("metSignificanceForCleaning", 3.0);
1053   iDesc.add<double>("metSignificanceForRejection", 4.0);
1054   iDesc.add<double>("metFactorForCleaning", 4.0);
1055   iDesc.add<double>("eventFractionForRejection", 0.8);
1056   iDesc.add<double>("metFactorForRejection", 4.0);
1057   iDesc.add<double>("metFactorForHighEta", 25.0);
1058   iDesc.add<double>("ptFactorForHighEta", 2.0);
1059   iDesc.add<double>("metFactorForFakes", 4.0);
1060   iDesc.add<double>("minMomentumForPunchThrough", 100.0);
1061   iDesc.add<double>("minEnergyForPunchThrough", 100.0);
1062   iDesc.add<double>("punchThroughFactor", 3.0);
1063   iDesc.add<double>("punchThroughMETFactor", 4.0);
1064   iDesc.add<double>("cosmicRejectionDistance", 1.0);
1065 }
[ 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!