HiJetAlgos/plugins/HiBadParticleCleaner.cc

0001 // system include files
0002 #include <memory>
0003
0004 // user include files
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 // class declaration
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   // ----------member data ---------------------------
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 // constructors and destructor
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 // member functions
0074 //
0075
0076 // ------------ method called on each new Event  ------------
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);  // mapping between the original PF and post-cleaning PF
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)  // muon cleaning
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)  //charged hadron cleaning
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   }  // end loop over pf candidates
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 //define this as a plug-in
0276 DEFINE_FWK_MODULE(HiBadParticleCleaner);