MuonAssociators/plugins/MuonCleanerBySegments.cc

0001
0002 //
0003 //
0004
0005 /**
0006   \class    modules::MuonCleanerBySegmentsT MuonCleanerBySegmentsT.h "MuonAnalysis/MuonAssociators/interface/MuonCleanerBySegmentsT.h"
0007   \brief    Removes duplicates from a muon collection using segment references.
0008
0009             This module removes duplicates from a muon collection using segment references.
0010
0011             All muons that don't pass the preselection are discarded first
0012
0013             Then, for each pair of muons that share at least a given fraction of segments,
0014             the worse one is flagged as ghost.
0015
0016             Finally, all muons that are not flagged as ghosts, or which pass a 'passthrough' selection,
0017             are saved in the output.
0018
0019   \author   Giovanni Petrucciani
0020   \version  $Id: MuonCleanerBySegments.cc,v 1.1 2012/08/11 13:00:50 gpetrucc Exp $
0021 */
0022
0023 #include "FWCore/Framework/interface/stream/EDProducer.h"
0024 #include "FWCore/Framework/interface/Event.h"
0025 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0026 #include "FWCore/Utilities/interface/InputTag.h"
0027
0028 #include "DataFormats/Math/interface/deltaR.h"
0029
0030 #include "DataFormats/MuonReco/interface/Muon.h"
0031 #include "DataFormats/MuonReco/interface/MuonSelectors.h"
0032 #include "DataFormats/PatCandidates/interface/Muon.h"
0033 #include "CommonTools/Utils/interface/StringCutObjectSelector.h"
0034
0035 namespace modules {
0036
0037   template <typename T>
0038   class MuonCleanerBySegmentsT : public edm::stream::EDProducer<> {
0039   public:
0040     explicit MuonCleanerBySegmentsT(const edm::ParameterSet &iConfig);
0041     ~MuonCleanerBySegmentsT() override {}
0042
0043     void produce(edm::Event &iEvent, const edm::EventSetup &iSetup) override;
0044
0045     bool isSameMuon(const T &mu1, const T &mu2) const {
0046       return (&mu1 == &mu2) ||
0047              (mu1.reco::Muon::innerTrack().isNonnull() ? mu1.reco::Muon::innerTrack() == mu2.reco::Muon::innerTrack()
0048                                                        : mu1.reco::Muon::outerTrack() == mu2.reco::Muon::outerTrack());
0049     }
0050     bool isBetterMuon(const T &mu1, const T &mu2) const;
0051
0052   private:
0053     /// Tokens for input collections
0054     edm::EDGetTokenT<edm::View<T> > srcToken_;
0055
0056     /// Preselection cut
0057     StringCutObjectSelector<T> preselection_;
0058     /// Always-accept cut
0059     StringCutObjectSelector<T> passthrough_;
0060
0061     /// Fraction of shared segments
0062     double sharedFraction_;
0063
0064     /// Use default criteria to choose the best muon
0065     bool defaultBestMuon_;
0066     /// Cut on the pair of objects together
0067     typedef std::pair<const reco::Muon *, const reco::Muon *> MuonPointerPair;
0068     StringCutObjectSelector<MuonPointerPair, true> bestMuonSelector_;
0069   };
0070
0071   template <>
0072   bool MuonCleanerBySegmentsT<pat::Muon>::isSameMuon(const pat::Muon &mu1, const pat::Muon &mu2) const;
0073 }  // namespace modules
0074
0075 template <typename T>
0076 modules::MuonCleanerBySegmentsT<T>::MuonCleanerBySegmentsT(const edm::ParameterSet &iConfig)
0077     : srcToken_(consumes<edm::View<T> >(iConfig.getParameter<edm::InputTag>("src"))),
0078       preselection_(iConfig.existsAs<std::string>("preselection") ? iConfig.getParameter<std::string>("preselection")
0079                                                                   : ""),
0080       passthrough_(iConfig.existsAs<std::string>("passthrough") ? iConfig.getParameter<std::string>("passthrough")
0081                                                                 : "0"),
0082       sharedFraction_(iConfig.getParameter<double>("fractionOfSharedSegments")),
0083       defaultBestMuon_(!iConfig.existsAs<std::string>("customArbitration")),
0084       bestMuonSelector_(defaultBestMuon_ ? std::string("") : iConfig.getParameter<std::string>("customArbitration")) {
0085   // this is the basic output (edm::Association is not generic)
0086   produces<std::vector<T> >();
0087 }
0088
0089 template <typename T>
0090 void modules::MuonCleanerBySegmentsT<T>::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) {
0091   using namespace edm;
0092   using namespace std;
0093
0094   Handle<View<T> > src;
0095   unique_ptr<vector<T> > out(new vector<T>());
0096
0097   iEvent.getByToken(srcToken_, src);
0098   unsigned int nsrc = src->size();
0099   out->reserve(nsrc);
0100   std::vector<int> good(nsrc, true);
0101   for (unsigned int i = 0; i < nsrc; ++i) {
0102     const T &mu1 = (*src)[i];
0103     if (!preselection_(mu1))
0104       good[i] = false;
0105     if (!good[i])
0106       continue;
0107     int nSegments1 = mu1.numberOfMatches(reco::Muon::SegmentArbitration);
0108     for (unsigned int j = i + 1; j < nsrc; ++j) {
0109       const T &mu2 = (*src)[j];
0110       if (isSameMuon(mu1, mu2))
0111         continue;
0112       if (!good[j] || !preselection_(mu2))
0113         continue;
0114       int nSegments2 = mu2.numberOfMatches(reco::Muon::SegmentArbitration);
0115       if (nSegments2 == 0 || nSegments1 == 0)
0116         continue;
0117       double sf = muon::sharedSegments(mu1, mu2) / std::min<double>(nSegments1, nSegments2);
0118       if (sf > sharedFraction_) {
0119         if (isBetterMuon(mu1, mu2)) {
0120           good[j] = false;
0121         } else {
0122           good[i] = false;
0123         }
0124       }
0125     }
0126   }
0127   for (unsigned int i = 0; i < nsrc; ++i) {
0128     const T &mu1 = (*src)[i];
0129     if (good[i] || passthrough_(mu1))
0130       out->push_back(mu1);
0131   }
0132   iEvent.put(std::move(out));
0133 }
0134
0135 template <typename T>
0136 bool modules::MuonCleanerBySegmentsT<T>::isBetterMuon(const T &mu1, const T &mu2) const {
0137   if (!defaultBestMuon_) {
0138     MuonPointerPair pair = {&mu1, &mu2};
0139     return bestMuonSelector_(pair);
0140   }
0141   if (mu2.track().isNull())
0142     return true;
0143   if (mu1.track().isNull())
0144     return false;
0145   if (mu1.isPFMuon() != mu2.isPFMuon())
0146     return mu1.isPFMuon();
0147   if (mu1.isGlobalMuon() != mu2.isGlobalMuon())
0148     return mu1.isGlobalMuon();
0149   if (mu1.charge() == mu2.charge() && deltaR2(mu1, mu2) < 0.0009) {
0150     return mu1.track()->ptError() / mu1.track()->pt() < mu2.track()->ptError() / mu2.track()->pt();
0151   } else {
0152     int nm1 = mu1.numberOfMatches(reco::Muon::SegmentArbitration);
0153     int nm2 = mu2.numberOfMatches(reco::Muon::SegmentArbitration);
0154     return (nm1 != nm2 ? nm1 > nm2 : mu1.pt() > mu2.pt());
0155   }
0156 }
0157
0158 template <>
0159 bool modules::MuonCleanerBySegmentsT<pat::Muon>::isSameMuon(const pat::Muon &mu1, const pat::Muon &mu2) const {
0160   return (&mu1 == &mu2) || (mu1.originalObjectRef() == mu2.originalObjectRef()) ||
0161          (mu1.reco::Muon::innerTrack().isNonnull() ? mu1.reco::Muon::innerTrack() == mu2.reco::Muon::innerTrack()
0162                                                    : mu1.reco::Muon::outerTrack() == mu2.reco::Muon::outerTrack());
0163 }
0164
0165 namespace modules {
0166   typedef modules::MuonCleanerBySegmentsT<reco::Muon> MuonCleanerBySegments;
0167   typedef modules::MuonCleanerBySegmentsT<pat::Muon> PATMuonCleanerBySegments;
0168 }  // namespace modules
0169
0170 #include "FWCore/Framework/interface/MakerMacros.h"
0171 using namespace modules;
0172 DEFINE_FWK_MODULE(MuonCleanerBySegments);
0173 DEFINE_FWK_MODULE(PATMuonCleanerBySegments);