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