File indexing completed on 2024-04-06 12:25:18
0001
0002 #include <memory>
0003
0004
0005 #include "FWCore/Framework/interface/Frameworkfwd.h"
0006 #include "FWCore/Framework/interface/global/EDProducer.h"
0007
0008 #include "FWCore/Framework/interface/Event.h"
0009 #include "FWCore/Framework/interface/MakerMacros.h"
0010
0011 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0012
0013 #include "DataFormats/MuonReco/interface/Muon.h"
0014 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0015 #include "DataFormats/MuonReco/interface/MuonSelectors.h"
0016 #include "DataFormats/TrackReco/interface/Track.h"
0017 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0018 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0019 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"
0020 #include "DataFormats/VertexReco/interface/Vertex.h"
0021 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0022
0023
0024
0025
0026 class HiBadParticleCleaner : public edm::global::EDProducer<> {
0027 public:
0028 explicit HiBadParticleCleaner(const edm::ParameterSet&);
0029 ~HiBadParticleCleaner() override = default;
0030
0031 private:
0032 void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0033
0034
0035 edm::EDGetTokenT<edm::View<reco::PFCandidate>> tokenPFCandidates_;
0036 edm::EDGetTokenT<reco::VertexCollection> tokenPV_;
0037
0038 const double minMuonPt_;
0039 const double minChargedHadronPt_;
0040 const double minMuonTrackRelPtErr_;
0041 const double maxSigLoose_;
0042 const double maxSigTight_;
0043 const double minCaloCompatibility_;
0044 const unsigned minTrackNHits_;
0045 const unsigned minPixelNHits_;
0046 const int minTrackerLayersForMuonLoose_;
0047 const int minTrackerLayersForMuonTight_;
0048 };
0049
0050
0051
0052
0053 HiBadParticleCleaner::HiBadParticleCleaner(const edm::ParameterSet& iConfig)
0054 : tokenPFCandidates_(consumes<edm::View<reco::PFCandidate>>(iConfig.getParameter<edm::InputTag>("PFCandidates"))),
0055 tokenPV_(consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("offlinePV"))),
0056 minMuonPt_(iConfig.getParameter<double>("minMuonPt")),
0057 minChargedHadronPt_(iConfig.getParameter<double>("minChargedHadronPt")),
0058 minMuonTrackRelPtErr_(iConfig.getParameter<double>("minMuonTrackRelPtErr")),
0059 maxSigLoose_(iConfig.getParameter<double>("maxSigLoose")),
0060 maxSigTight_(iConfig.getParameter<double>("maxSigTight")),
0061 minCaloCompatibility_(iConfig.getParameter<double>("minCaloCompatibility")),
0062 minTrackNHits_(iConfig.getParameter<uint>("minTrackNHits")),
0063 minPixelNHits_(iConfig.getParameter<uint>("minPixelNHits")),
0064 minTrackerLayersForMuonLoose_(iConfig.getParameter<int>("minTrackerLayersForMuonLoose")),
0065 minTrackerLayersForMuonTight_(iConfig.getParameter<int>("minTrackerLayersForMuonTight")) {
0066 produces<bool>();
0067 produces<reco::PFCandidateCollection>();
0068 produces<reco::PFCandidateCollection>("removed");
0069 produces<edm::ValueMap<reco::PFCandidateRef>>();
0070 }
0071
0072
0073
0074
0075
0076
0077 void HiBadParticleCleaner::produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup&) const {
0078 using namespace std;
0079 using namespace edm;
0080
0081 typedef View<reco::PFCandidate> CandidateView;
0082 Handle<CandidateView> pfCandidates;
0083 iEvent.getByToken(tokenPFCandidates_, pfCandidates);
0084
0085 const reco::VertexCollection* recoVertices;
0086 edm::Handle<reco::VertexCollection> vertexCollection;
0087 iEvent.getByToken(tokenPV_, vertexCollection);
0088 recoVertices = vertexCollection.product();
0089
0090 auto pOutputCandidateCollection = std::make_unique<reco::PFCandidateCollection>();
0091 auto pBadCandidateCollection = std::make_unique<reco::PFCandidateCollection>();
0092
0093 bool foundBadCandidate = false;
0094
0095 size_t n = pfCandidates->size();
0096 std::vector<int> candidateIndexMapper(n, 0);
0097 size_t iPF;
0098 for (iPF = 0; iPF < n; iPF++) {
0099 const reco::PFCandidate& pfCandidate = pfCandidates->at(iPF);
0100 if (pfCandidate.particleId() == reco::PFCandidate::ParticleType::mu)
0101 {
0102 if (pfCandidate.pt() > minMuonPt_) {
0103 if (!pfCandidate.muonRef()->isGlobalMuon() || !pfCandidate.muonRef()->isTrackerMuon() ||
0104 !pfCandidate.trackRef().isNonnull()) {
0105 foundBadCandidate = true;
0106 continue;
0107 }
0108 reco::TrackRef track = pfCandidate.trackRef();
0109
0110 if (track->ptError() / track->pt() > minMuonTrackRelPtErr_ || track->pt() < pfCandidate.pt() / 2.) {
0111 foundBadCandidate = true;
0112 continue;
0113 }
0114
0115 if (track->algo() == reco::TrackBase::muonSeededStepInOut ||
0116 track->algo() == reco::TrackBase::muonSeededStepOutIn ||
0117 track->originalAlgo() == reco::TrackBase::muonSeededStepInOut ||
0118 track->originalAlgo() == reco::TrackBase::muonSeededStepOutIn ||
0119 track->hitPattern().trackerLayersWithMeasurement() < minTrackerLayersForMuonLoose_) {
0120 const reco::Vertex& vtx = (*recoVertices)[0];
0121 float bestVzError = vtx.zError();
0122 const math::XYZPoint& bestVtx(vtx.position());
0123 math::Error<3>::type vtx_cov = vtx.covariance();
0124 float dz = std::abs(track->dz(bestVtx));
0125 float dxy = std::abs(track->dxy(bestVtx));
0126 float dzError2 = track->dzError() * track->dzError() + bestVzError * bestVzError;
0127 float dxyError = track->dxyError(bestVtx, vtx_cov);
0128
0129 float dzSig2 = dz * dz / dzError2;
0130 float dxySig2 = dxy * dxy / dxyError / dxyError;
0131
0132 float sig3d = sqrt(dzSig2 + dxySig2);
0133
0134 if (sig3d > maxSigLoose_) {
0135 pBadCandidateCollection->push_back(pfCandidate);
0136 candidateIndexMapper[iPF] = -1 * (pBadCandidateCollection->size());
0137 foundBadCandidate = true;
0138 continue;
0139 }
0140
0141 if (track->pt() < pfCandidate.pt() / 1.5 || track->pt() > pfCandidate.pt() * 1.5) {
0142 foundBadCandidate = true;
0143 pBadCandidateCollection->push_back(pfCandidate);
0144 candidateIndexMapper[iPF] = -1 * (pBadCandidateCollection->size());
0145 continue;
0146 }
0147 if (track->originalAlgo() == reco::TrackBase::muonSeededStepOutIn &&
0148 track->hitPattern().trackerLayersWithMeasurement() < minTrackerLayersForMuonTight_) {
0149 foundBadCandidate = true;
0150 pBadCandidateCollection->push_back(pfCandidate);
0151 candidateIndexMapper[iPF] = -1 * (pBadCandidateCollection->size());
0152 continue;
0153 }
0154 }
0155 }
0156 } else if (pfCandidate.particleId() == reco::PFCandidate::ParticleType::h)
0157 {
0158 if (pfCandidate.pt() > minChargedHadronPt_) {
0159 reco::TrackRef track = pfCandidate.trackRef();
0160
0161 unsigned nHits = track->numberOfValidHits();
0162 unsigned nPixelHits = track->hitPattern().numberOfValidPixelHits();
0163
0164 if ((nHits < minTrackNHits_ && nPixelHits < minPixelNHits_) || nHits == 3) {
0165 foundBadCandidate = true;
0166 pBadCandidateCollection->push_back(pfCandidate);
0167 candidateIndexMapper[iPF] = -1 * (pBadCandidateCollection->size());
0168 continue;
0169 }
0170
0171 const reco::Vertex& vtx = (*recoVertices)[0];
0172 float bestVzError = vtx.zError();
0173 const math::XYZPoint& bestVtx(vtx.position());
0174 math::Error<3>::type vtx_cov = vtx.covariance();
0175 float dz = std::abs(track->dz(bestVtx));
0176 float dxy = std::abs(track->dxy(bestVtx));
0177 float dzError2 = track->dzError() * track->dzError() + bestVzError * bestVzError;
0178 float dxyError = track->dxyError(bestVtx, vtx_cov);
0179
0180 float dzSig2 = dz * dz / dzError2;
0181 float dxySig2 = dxy * dxy / dxyError / dxyError;
0182
0183 float sig3d = sqrt(dzSig2 + dxySig2);
0184
0185 if (sig3d > maxSigLoose_) {
0186 foundBadCandidate = true;
0187 pBadCandidateCollection->push_back(pfCandidate);
0188 candidateIndexMapper[iPF] = -1 * (pBadCandidateCollection->size());
0189 continue;
0190 }
0191
0192 if (sig3d > maxSigTight_ && nHits < minTrackNHits_) {
0193 foundBadCandidate = true;
0194 pBadCandidateCollection->push_back(pfCandidate);
0195 candidateIndexMapper[iPF] = -1 * (pBadCandidateCollection->size());
0196 continue;
0197 }
0198
0199 if (track->algo() == reco::TrackBase::muonSeededStepInOut ||
0200 track->algo() == reco::TrackBase::muonSeededStepOutIn ||
0201 track->originalAlgo() == reco::TrackBase::muonSeededStepInOut ||
0202 track->originalAlgo() == reco::TrackBase::muonSeededStepOutIn) {
0203 if (sig3d > maxSigLoose_) {
0204 foundBadCandidate = true;
0205 pBadCandidateCollection->push_back(pfCandidate);
0206 candidateIndexMapper[iPF] = -1 * (pBadCandidateCollection->size());
0207 continue;
0208 }
0209
0210 if (nHits < minTrackNHits_) {
0211 foundBadCandidate = true;
0212 pBadCandidateCollection->push_back(pfCandidate);
0213 candidateIndexMapper[iPF] = -1 * (pBadCandidateCollection->size());
0214 continue;
0215 }
0216 }
0217
0218 double caloEnergy = pfCandidate.ecalEnergy() + pfCandidate.hcalEnergy();
0219
0220 if (caloEnergy < track->p() * minCaloCompatibility_) {
0221 if (sig3d > maxSigTight_) {
0222 foundBadCandidate = true;
0223 pBadCandidateCollection->push_back(pfCandidate);
0224 candidateIndexMapper[iPF] = -1 * (pBadCandidateCollection->size());
0225 continue;
0226 }
0227
0228 if (nHits < minTrackNHits_) {
0229 foundBadCandidate = true;
0230 pBadCandidateCollection->push_back(pfCandidate);
0231 candidateIndexMapper[iPF] = -1 * (pBadCandidateCollection->size());
0232 continue;
0233 }
0234
0235 if (nPixelHits < minPixelNHits_) {
0236 foundBadCandidate = true;
0237 pBadCandidateCollection->push_back(pfCandidate);
0238 candidateIndexMapper[iPF] = -1 * (pBadCandidateCollection->size());
0239 continue;
0240 }
0241 }
0242 }
0243 }
0244
0245 pOutputCandidateCollection->push_back(pfCandidate);
0246 candidateIndexMapper[iPF] = (pOutputCandidateCollection->size());
0247 }
0248
0249 bool pass = !foundBadCandidate;
0250
0251 edm::OrphanHandle<std::vector<reco::PFCandidate>> newpf = iEvent.put(std::move(pOutputCandidateCollection));
0252 edm::OrphanHandle<std::vector<reco::PFCandidate>> badpf = iEvent.put(std::move(pBadCandidateCollection), "removed");
0253
0254 iEvent.put(std::make_unique<bool>(pass));
0255
0256 std::unique_ptr<edm::ValueMap<reco::PFCandidateRef>> pf2pf(new edm::ValueMap<reco::PFCandidateRef>());
0257 edm::ValueMap<reco::PFCandidateRef>::Filler filler(*pf2pf);
0258
0259 std::vector<reco::PFCandidateRef> refs;
0260 refs.reserve(n);
0261
0262 for (iPF = 0; iPF < n; ++iPF) {
0263 if (candidateIndexMapper[iPF] > 0) {
0264 refs.push_back(reco::PFCandidateRef(newpf, candidateIndexMapper[iPF] - 1));
0265 } else if (candidateIndexMapper[iPF] < 0) {
0266 refs.push_back(reco::PFCandidateRef(badpf, -candidateIndexMapper[iPF] - 1));
0267 }
0268 }
0269 filler.insert(pfCandidates, refs.begin(), refs.end());
0270
0271 filler.fill();
0272 iEvent.put(std::move(pf2pf));
0273 }
0274
0275
0276 DEFINE_FWK_MODULE(HiBadParticleCleaner);