MuonAssociators/plugins/MuonFSRProducer.cc

0001 // system include files
0002 #include <memory>
0003
0004 // CMSSW include files
0005 #include "DataFormats/Candidate/interface/Candidate.h"
0006 #include "DataFormats/Math/interface/LorentzVector.h"
0007 #include "DataFormats/PatCandidates/interface/Electron.h"
0008 #include "DataFormats/PatCandidates/interface/GenericParticle.h"
0009 #include "DataFormats/PatCandidates/interface/Muon.h"
0010 #include "DataFormats/PatCandidates/interface/PackedCandidate.h"
0011 #include "FWCore/Framework/interface/Event.h"
0012 #include "FWCore/Framework/interface/Frameworkfwd.h"
0013 #include "FWCore/Framework/interface/MakerMacros.h"
0014 #include "FWCore/Framework/interface/global/EDProducer.h"
0015 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0016 #include "FWCore/Utilities/interface/StreamID.h"
0017
0018 //
0019 // class declaration
0020 //
0021
0022 class MuonFSRProducer : public edm::global::EDProducer<> {
0023 public:
0024   explicit MuonFSRProducer(const edm::ParameterSet& iConfig)
0025       :
0026
0027         pfcands_{consumes<pat::PackedCandidateCollection>(iConfig.getParameter<edm::InputTag>("packedPFCandidates"))},
0028         electrons_{consumes<pat::ElectronCollection>(iConfig.getParameter<edm::InputTag>("slimmedElectrons"))},
0029         muons_{consumes<edm::View<reco::Muon>>(iConfig.getParameter<edm::InputTag>("muons"))},
0030         ptCut(iConfig.getParameter<double>("muonPtMin")),
0031         etaCut(iConfig.getParameter<double>("muonEtaMax")),
0032         photonPtCut(iConfig.getParameter<double>("photonPtMin")),
0033         drEtCut(iConfig.getParameter<double>("deltaROverEt2Max")),
0034         isoCut(iConfig.getParameter<double>("isolation")) {
0035     produces<std::vector<pat::GenericParticle>>();
0036   }
0037   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0038     edm::ParameterSetDescription desc;
0039     desc.add<edm::InputTag>("packedPFCandidates", edm::InputTag("packedPFCandidates"))
0040         ->setComment("packed pf candidates where to look for photons");
0041     desc.add<edm::InputTag>("slimmedElectrons", edm::InputTag("slimmedElectrons"))
0042         ->setComment(
0043             "electrons to check for footprint, the electron collection must have proper linking with the "
0044             "packedCandidate collection");
0045     desc.add<edm::InputTag>("muons", edm::InputTag("slimmedMuons"))
0046         ->setComment("collection of muons to correct for FSR ");
0047     desc.add<double>("muonPtMin", 20.)->setComment("minimum pt of the muon to look for a near photon");
0048     desc.add<double>("muonEtaMax", 2.4)->setComment("max eta of the muon to look for a near photon");
0049     desc.add<double>("photonPtMin", 2.0)->setComment("minimum photon Pt");
0050     desc.add<double>("deltaROverEt2Max", 0.05)->setComment("max ratio of deltsR(mu,photon) over et2 of the photon");
0051     desc.add<double>("isolation", 2.0)->setComment("relative isolation cut");
0052
0053     descriptions.addWithDefaultLabel(desc);
0054   }
0055   ~MuonFSRProducer() override {}
0056
0057 private:
0058   void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0059
0060   double computeRelativeIsolation(const pat::PackedCandidate& photon,
0061                                   const pat::PackedCandidateCollection& pfcands,
0062                                   const double& isoConeMax,
0063                                   const double& isoConeMin) const;
0064
0065   // ----------member data ---------------------------
0066   const edm::EDGetTokenT<pat::PackedCandidateCollection> pfcands_;
0067   const edm::EDGetTokenT<pat::ElectronCollection> electrons_;
0068   const edm::EDGetTokenT<edm::View<reco::Muon>> muons_;
0069   float ptCut;
0070   float etaCut;
0071   float photonPtCut;
0072   float drEtCut;
0073   float isoCut;
0074 };
0075
0076 void MuonFSRProducer::produce(edm::StreamID streamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0077   using namespace std;
0078
0079   edm::Handle<pat::PackedCandidateCollection> pfcands;
0080   iEvent.getByToken(pfcands_, pfcands);
0081   edm::Handle<edm::View<reco::Muon>> muons;
0082   iEvent.getByToken(muons_, muons);
0083   edm::Handle<pat::ElectronCollection> electrons;
0084   iEvent.getByToken(electrons_, electrons);
0085
0086   auto fsrPhotons = std::make_unique<std::vector<pat::GenericParticle>>();
0087   // loop over all muons
0088   for (auto muon = muons->begin(); muon != muons->end(); ++muon) {
0089     int photonPosition = -1;
0090     double distance_metric_min = -1;
0091     // minimum muon pT
0092     if (muon->pt() < ptCut)
0093       continue;
0094     // maximum muon eta
0095     if (abs(muon->eta()) > etaCut)
0096       continue;
0097
0098     // for each muon, loop over all pf cadidates
0099     for (auto iter_pf = pfcands->begin(); iter_pf != pfcands->end(); iter_pf++) {
0100       auto const& pc = *iter_pf;
0101
0102       // consider only photons
0103       if (abs(pc.pdgId()) != 22)
0104         continue;
0105       // minimum pT cut
0106       if (pc.pt() < photonPtCut)
0107         continue;
0108
0109       // eta requirements
0110       if (abs(pc.eta()) > 1.4 and (abs(pc.eta()) < 1.6))
0111         continue;
0112       if (abs(pc.eta()) > 2.5)
0113         continue;
0114
0115       // 0.0001 < DeltaR(photon,muon) < 0.5 requirement
0116       double dRPhoMu = deltaR(muon->eta(), muon->phi(), pc.eta(), pc.phi());
0117       if (dRPhoMu < 0.0001)
0118         continue;
0119       if (dRPhoMu > 0.5)
0120         continue;
0121
0122       bool skipPhoton = false;
0123       bool closest = true;
0124
0125       for (auto othermuon = muons->begin(); othermuon != muons->end(); ++othermuon) {
0126         if (othermuon->pt() < ptCut or abs(othermuon->eta()) > etaCut or muon == othermuon)
0127           continue;
0128         double dRPhoMuOther = deltaR(othermuon->eta(), othermuon->phi(), pc.eta(), pc.phi());
0129         if (dRPhoMuOther < dRPhoMu) {
0130           closest = false;
0131           break;
0132         }
0133       }
0134       if (!closest)
0135         continue;
0136
0137       // Check that is not in footprint of an electron
0138       pat::PackedCandidateRef pfcandRef = pat::PackedCandidateRef(pfcands, iter_pf - pfcands->begin());
0139
0140       for (auto electrons_iter = electrons->begin(); electrons_iter != electrons->end(); ++electrons_iter) {
0141         for (auto const& cand : electrons_iter->associatedPackedPFCandidates()) {
0142           if (!cand.isAvailable())
0143             continue;
0144           if (cand.id() != pfcandRef.id())
0145             throw cms::Exception("Configuration")
0146                 << "The electron associatedPackedPFCandidates item does not have "
0147                 << "the same ID of packed candidate collection used for cleaning the electron footprint: " << cand.id()
0148                 << " (" << pfcandRef.id() << ")\n";
0149           if (cand.key() == pfcandRef.key()) {
0150             skipPhoton = true;
0151             break;
0152           }
0153         }
0154         if (skipPhoton)
0155           break;
0156       }
0157
0158       if (skipPhoton)
0159         continue;
0160
0161       // use only isolated photons (very loose prelection can be tightened on analysis level)
0162       float photon_relIso03 = computeRelativeIsolation(pc, *pfcands, 0.3, 0.0001);
0163       if (photon_relIso03 > isoCut)
0164         continue;
0165       double metric = deltaR(muon->eta(), muon->phi(), pc.eta(), pc.phi()) / (pc.pt() * pc.pt());
0166       if (metric > drEtCut)
0167         continue;
0168       fsrPhotons->push_back(pat::GenericParticle(pc));
0169       fsrPhotons->back().addUserFloat("relIso03", photon_relIso03);  // isolation, no CHS
0170       fsrPhotons->back().addUserCand("associatedMuon", reco::CandidatePtr(muons, muon - muons->begin()));
0171       fsrPhotons->back().addUserFloat("dROverEt2", metric);  // dR/et2 to the closest muon
0172
0173       // FSR photon defined as the one with minimum value of DeltaR/Et^2
0174       if (photonPosition == -1 or metric < distance_metric_min) {
0175         distance_metric_min = metric;
0176         photonPosition = fsrPhotons->size() - 1;
0177       }
0178     }
0179   }
0180
0181   iEvent.put(std::move(fsrPhotons));
0182 }
0183
0184 double MuonFSRProducer::computeRelativeIsolation(const pat::PackedCandidate& photon,
0185                                                  const pat::PackedCandidateCollection& pfcands,
0186                                                  const double& isoConeMax,
0187                                                  const double& isoConeMin) const {
0188   double ptsum = 0;
0189
0190   for (const auto& pfcand : pfcands) {
0191     // Isolation cone requirement
0192     double dRIsoCone = deltaR(photon.eta(), photon.phi(), pfcand.eta(), pfcand.phi());
0193     if (dRIsoCone > isoConeMax)
0194       continue;
0195     if (dRIsoCone < isoConeMin)
0196       continue;
0197
0198     if (pfcand.charge() != 0 && abs(pfcand.pdgId()) == 211 && pfcand.pt() > 0.2) {
0199       if (dRIsoCone > 0.0001)
0200         ptsum += pfcand.pt();
0201     } else if (pfcand.charge() == 0 && (abs(pfcand.pdgId()) == 22 || abs(pfcand.pdgId()) == 130) && pfcand.pt() > 0.5) {
0202       if (dRIsoCone > 0.01)
0203         ptsum += pfcand.pt();
0204     }
0205   }
0206
0207   return ptsum / photon.pt();
0208 }
0209
0210 //define this as a plug-in
0211 DEFINE_FWK_MODULE(MuonFSRProducer);