Merge pull request #36323 from vpicco/NANO_FSRForEle

Include FSR for electrons (in addition to muons)

CommonTools/RecoUtils/plugins/LeptonFSRProducer.cc

@@ -0,0 +1,286 @@
+/** \class LeptonFSRProducer
+ * Search for FSR photons for muons and electrons.
+ *
+ * Photon candidates are searched among the "packedPFCandidates" collection with the specified cuts, and are required to be isolated 
+ * (relIso, with a cone of 0.3) and not to be in the footprint of all electrons in the "electrons" collection.
+ * Each photon is matched by DeltaR to the closest among all muons and electrons and stored if passing dR/ET^2<deltaROverEt2Max.
+ * In addition ValueMaps are stored, with links to one photon per muon/electron. For this purpose, if more than a photon
+ * is matched to a lepton, the lowest-DR/ET^2 is chosen.
+ *
+ */
+
+#include <memory>
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+#include "FWCore/Framework/interface/global/EDProducer.h"
+
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/Utilities/interface/StreamID.h"
+
+#include "DataFormats/Candidate/interface/Candidate.h"
+#include "DataFormats/PatCandidates/interface/PackedCandidate.h"
+#include "DataFormats/PatCandidates/interface/GenericParticle.h"
+#include "DataFormats/Math/interface/LorentzVector.h"
+
+#include "DataFormats/PatCandidates/interface/Muon.h"
+#include "DataFormats/PatCandidates/interface/Electron.h"
+#include "DataFormats/Common/interface/ValueMap.h"
+
+class LeptonFSRProducer : public edm::global::EDProducer<> {
+public:
+  explicit LeptonFSRProducer(const edm::ParameterSet& iConfig)
+      : pfcands_{consumes<pat::PackedCandidateCollection>(iConfig.getParameter<edm::InputTag>("packedPFCandidates"))},
+        electronsForVeto_{consumes<pat::ElectronCollection>(iConfig.getParameter<edm::InputTag>("slimmedElectrons"))},
+        muons_{consumes<edm::View<reco::Muon>>(iConfig.getParameter<edm::InputTag>("muons"))},
+        electrons_{consumes<edm::View<reco::GsfElectron>>(iConfig.getParameter<edm::InputTag>("electrons"))},
+        ptCutMu(iConfig.getParameter<double>("muonPtMin")),
+        etaCutMu(iConfig.getParameter<double>("muonEtaMax")),
+        ptCutE(iConfig.getParameter<double>("elePtMin")),
+        etaCutE(iConfig.getParameter<double>("eleEtaMax")),
+        photonPtCut(iConfig.getParameter<double>("photonPtMin")),
+        drEtCut(iConfig.getParameter<double>("deltaROverEt2Max")),
+        isoCut(iConfig.getParameter<double>("isolation")),
+        drSafe(0.0001) {
+    produces<std::vector<pat::GenericParticle>>();
+    produces<edm::ValueMap<int>>("muFsrIndex");
+    produces<edm::ValueMap<int>>("eleFsrIndex");
+  }
+  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+    edm::ParameterSetDescription desc;
+    desc.add<edm::InputTag>("packedPFCandidates", edm::InputTag("packedPFCandidates"))
+        ->setComment("packed pf candidates where to look for photons");
+    desc.add<edm::InputTag>("slimmedElectrons", edm::InputTag("slimmedElectrons"))
+        ->setComment(
+            "electrons to check for footprint, the electron collection must have proper linking with the "
+            "packedCandidate collection");
+    desc.add<edm::InputTag>("muons", edm::InputTag("slimmedMuons"))
+        ->setComment("collection of muons to match with FSR ");
+    desc.add<edm::InputTag>("electrons", edm::InputTag("slimmedElectrons"))
+        ->setComment("collection of electrons to match with FSR ");
+    desc.add<double>("muonPtMin", 3.)->setComment("minimum pt of the muon to look for a near photon");
+    desc.add<double>("muonEtaMax", 2.4)->setComment("max eta of the muon to look for a near photon");
+    desc.add<double>("elePtMin", 5.)->setComment("minimum pt of the electron to look for a near photon");
+    desc.add<double>("eleEtaMax", 2.5)->setComment("max eta of the electron to look for a near photon");
+    desc.add<double>("photonPtMin", 2.0)->setComment("minimum photon Pt");
+    desc.add<double>("deltaROverEt2Max", 0.05)->setComment("max ratio of deltaR(lep,photon) over et2 of the photon");
+    desc.add<double>("isolation", 2.0)->setComment("photon relative isolation cut");
+
+    descriptions.addWithDefaultLabel(desc);
+  }
+  ~LeptonFSRProducer() override = default;
+
+private:
+  void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
+
+  double computeRelativeIsolation(const pat::PackedCandidate& photon,
+                                  const pat::PackedCandidateCollection& pfcands,
+                                  const double& isoConeMax2,
+                                  const double& isoConeMin2) const;
+
+  bool electronFootprintVeto(pat::PackedCandidateRef& pfcandRef,
+                             edm::Handle<pat::ElectronCollection> electronsForVeto) const;
+
+  // ----------member data ---------------------------
+  const edm::EDGetTokenT<pat::PackedCandidateCollection> pfcands_;
+  const edm::EDGetTokenT<pat::ElectronCollection> electronsForVeto_;
+  const edm::EDGetTokenT<edm::View<reco::Muon>> muons_;
+  const edm::EDGetTokenT<edm::View<reco::GsfElectron>> electrons_;
+  const double ptCutMu;
+  const double etaCutMu;
+  const double ptCutE;
+  const double etaCutE;
+  const double photonPtCut;
+  const double drEtCut;
+  const double isoCut;
+  const double drSafe;
+};
+
+void LeptonFSRProducer::produce(edm::StreamID streamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
+  using namespace std;
+
+  edm::Handle<pat::PackedCandidateCollection> pfcands;
+  iEvent.getByToken(pfcands_, pfcands);
+  edm::Handle<edm::View<reco::Muon>> muons;
+  iEvent.getByToken(muons_, muons);
+  edm::Handle<edm::View<reco::GsfElectron>> electrons;
+  iEvent.getByToken(electrons_, electrons);
+  edm::Handle<pat::ElectronCollection> electronsForVeto;
+  iEvent.getByToken(electronsForVeto_, electronsForVeto);
+
+  // The output collection of FSR photons
+  auto fsrPhotons = std::make_unique<std::vector<pat::GenericParticle>>();
+
+  std::vector<int> muPhotonIdxs(muons->size(), -1);
+  std::vector<double> muPhotonDRET2(muons->size(), 1e9);
+
+  std::vector<int> elePhotonIdxs(electrons->size(), -1);
+  std::vector<double> elePhotonDRET2(electrons->size(), 1e9);
+
+  //----------------------
+  // Loop on photon candidates
+  //----------------------
+
+  for (auto pc = pfcands->begin(); pc != pfcands->end(); pc++) {
+    // consider only photons, with pT and eta cuts
+    if (abs(pc->pdgId()) != 22 || pc->pt() < photonPtCut || abs(pc->eta()) > 2.5)
+      continue;
+
+    //------------------------------------------------------
+    // Get the closest lepton
+    //------------------------------------------------------
+    double dRMin(0.5);
+    int closestMu = -1;
+    int closestEle = -1;
+    double photon_relIso03 = 1e9;  // computed only if necessary
+    bool skipPhoton = false;
+
+    for (auto muon = muons->begin(); muon != muons->end(); ++muon) {
+      if (muon->pt() < ptCutMu || std::abs(muon->eta()) > etaCutMu)
+        continue;
+
+      int muonIdx = muon - muons->begin();
+      double dR = deltaR(muon->eta(), muon->phi(), pc->eta(), pc->phi());
+      if (dR < dRMin && dR > drSafe && dR < drEtCut * pc->pt() * pc->pt()) {
+        // Check if photon is isolated
+        photon_relIso03 = computeRelativeIsolation(*pc, *pfcands, 0.3 * 0.3, drSafe * drSafe);
+        if (photon_relIso03 > isoCut) {
+          skipPhoton = true;
+          break;  // break loop on muons -> photon will be skipped
+        }
+        // Check that photon is not in footprint of an electron
+        pat::PackedCandidateRef pfcandRef = pat::PackedCandidateRef(pfcands, pc - pfcands->begin());
+        skipPhoton = electronFootprintVeto(pfcandRef, electronsForVeto);
+        if (skipPhoton)
+          break;  // break loop on muons -> photon will be skipped
+
+        // Candidate matching
+        dRMin = dR;
+        closestMu = muonIdx;
+      }
+    }  // end of loop on muons
+
+    if (skipPhoton)
+      continue;  // photon does not pass iso or ele footprint veto; do not look for electrons
+
+    for (auto ele = electrons->begin(); ele != electrons->end(); ++ele) {
+      if (ele->pt() < ptCutE || std::abs(ele->eta()) > etaCutE)
+        continue;
+
+      int eleIdx = ele - electrons->begin();
+      double dR = deltaR(ele->eta(), ele->phi(), pc->eta(), pc->phi());
+      if (dR < dRMin && dR > drSafe && dR < drEtCut * pc->pt() * pc->pt()) {
+        // Check if photon is isolated (no need to recompute iso if already done for muons above)
+        if (photon_relIso03 > 1e8) {
+          photon_relIso03 = computeRelativeIsolation(*pc, *pfcands, 0.3 * 0.3, drSafe * drSafe);
+        }
+        if (photon_relIso03 > isoCut) {
+          break;  // break loop on electrons -> photon will be skipped
+        }
+        // Check that photon is not in footprint of an electron
+        pat::PackedCandidateRef pfcandRef = pat::PackedCandidateRef(pfcands, pc - pfcands->begin());
+        if (electronFootprintVeto(pfcandRef, electronsForVeto)) {
+          break;  // break loop on electrons -> photon will be skipped
+        }
+
+        // Candidate matching
+        dRMin = dR;
+        closestEle = eleIdx;
+        closestMu = -1;  // reset match to muons
+      }
+    }  // end loop on electrons
+
+    if (closestMu >= 0 || closestEle >= 0) {
+      // Add FSR photon to the output collection
+      double dRET2 = dRMin / pc->pt() / pc->pt();
+      int iPhoton = fsrPhotons->size();
+      fsrPhotons->push_back(pat::GenericParticle(*pc));
+      fsrPhotons->back().addUserFloat("relIso03", photon_relIso03);
+      fsrPhotons->back().addUserFloat("dROverEt2", dRET2);
+
+      if (closestMu >= 0) {
+        fsrPhotons->back().addUserCand("associatedMuon", reco::CandidatePtr(muons, closestMu));
+        // Store the backlink to the photon: choose the lowest-dRET2 photon for each mu...
+        if (dRET2 < muPhotonDRET2[closestMu]) {
+          muPhotonDRET2[closestMu] = dRET2;
+          muPhotonIdxs[closestMu] = iPhoton;
+        }
+      } else if (closestEle >= 0) {
+        // ...and same for eles
+        fsrPhotons->back().addUserCand("associatedElectron", reco::CandidatePtr(electrons, closestEle));
+        if (dRET2 < elePhotonDRET2[closestEle]) {
+          elePhotonDRET2[closestEle] = dRET2;
+          elePhotonIdxs[closestEle] = iPhoton;
+        }
+      }
+    }
+  }  // end of loop over pfCands
+
+  edm::OrphanHandle<std::vector<pat::GenericParticle>> oh = iEvent.put(std::move(fsrPhotons));
+
+  {
+    std::unique_ptr<edm::ValueMap<int>> bareIdx(new edm::ValueMap<int>());
+    edm::ValueMap<int>::Filler fillerBareIdx(*bareIdx);
+    fillerBareIdx.insert(muons, muPhotonIdxs.begin(), muPhotonIdxs.end());
+    fillerBareIdx.fill();
+    iEvent.put(std::move(bareIdx), "muFsrIndex");
+  }
+
+  {
+    std::unique_ptr<edm::ValueMap<int>> bareIdx(new edm::ValueMap<int>());
+    edm::ValueMap<int>::Filler fillerBareIdx(*bareIdx);
+    fillerBareIdx.insert(electrons, elePhotonIdxs.begin(), elePhotonIdxs.end());
+    fillerBareIdx.fill();
+    iEvent.put(std::move(bareIdx), "eleFsrIndex");
+  }
+}
+
+double LeptonFSRProducer::computeRelativeIsolation(const pat::PackedCandidate& photon,
+                                                   const pat::PackedCandidateCollection& pfcands,
+                                                   const double& isoConeMax2,
+                                                   const double& isoConeMin2) const {
+  double ptsum = 0;
+
+  for (const auto& pfcand : pfcands) {
+    // Isolation cone
+    double dRIsoCone2 = deltaR2(photon.eta(), photon.phi(), pfcand.eta(), pfcand.phi());
+    if (dRIsoCone2 > isoConeMax2 || dRIsoCone2 < isoConeMin2)
+      continue;
+
+    // Charged hadrons
+    if (pfcand.charge() != 0 && abs(pfcand.pdgId()) == 211 && pfcand.pt() > 0.2 && dRIsoCone2 > drSafe * drSafe) {
+      ptsum += pfcand.pt();
+      // Neutral hadrons + photons
+    } else if (pfcand.charge() == 0 && (abs(pfcand.pdgId()) == 22 || abs(pfcand.pdgId()) == 130) && pfcand.pt() > 0.5 &&
+               dRIsoCone2 > 0.01 * 0.01) {
+      ptsum += pfcand.pt();
+    }
+  }
+
+  return ptsum / photon.pt();
+}
+
+bool LeptonFSRProducer::electronFootprintVeto(pat::PackedCandidateRef& pfcandRef,
+                                              edm::Handle<pat::ElectronCollection> electronsForVeto) const {
+  bool skipPhoton = false;
+  for (auto electrons_iter = electronsForVeto->begin(); electrons_iter != electronsForVeto->end(); ++electrons_iter) {
+    for (auto const& cand : electrons_iter->associatedPackedPFCandidates()) {
+      if (!cand.isAvailable())
+        continue;
+      if (cand.id() != pfcandRef.id())
+        throw cms::Exception("Configuration")
+            << "The electron associatedPackedPFCandidates item does not have "
+            << "the same ID of packed candidate collection used for cleaning the electron footprint: " << cand.id()
+            << " (" << pfcandRef.id() << ")\n";
+      if (cand.key() == pfcandRef.key()) {
+        skipPhoton = true;
+        break;
+      }
+    }
+  }
+  return skipPhoton;
+}
+
+//define this as a plug-in
+DEFINE_FWK_MODULE(LeptonFSRProducer);

PhysicsTools/NanoAOD/python/electrons_cff.py

@@ -395,6 +395,7 @@ def _get_bitmapVIDForEle_docstring(modules,WorkingPoints):
     ),
     externalVariables = cms.PSet(
         mvaTTH = ExtVar(cms.InputTag("electronMVATTH"),float, doc="TTH MVA lepton ID score",precision=14),
+        fsrPhotonIdx = ExtVar(cms.InputTag("leptonFSRphotons:eleFsrIndex"),int, doc="Index of the lowest-dR/ET2 among associated FSR photons"),
     ),
 )
 
@@ -554,7 +555,8 @@ def _get_bitmapVIDForEle_docstring(modules,WorkingPoints):
 )
 
 (run3_nanoAOD_devel).toReplaceWith(electronTablesTask, electronTablesTask.copyAndExclude([electronMVATTH]))
-(run3_nanoAOD_devel).toModify(electronTable, externalVariables = None)
+(run3_nanoAOD_devel).toModify(electronTable, externalVariables = cms.PSet(fsrPhotonIdx = ExtVar(cms.InputTag("leptonFSRphotons:eleFsrIndex"),int, doc="Index of the  lowest-dR/ET2 among associated FSR photons")),
+)
 ##### end TEMPORARY Run3
 
 

PhysicsTools/NanoAOD/python/fsrPhotons_cff.py

@@ -0,0 +1,27 @@
+import FWCore.ParameterSet.Config as cms
+from PhysicsTools.NanoAOD.common_cff import *
+
+from CommonTools.RecoUtils.leptonFSRProducer_cfi import leptonFSRProducer
+leptonFSRphotons = leptonFSRProducer.clone(
+  packedPFCandidates = "packedPFCandidates",
+  slimmedElectrons = "slimmedElectrons", #for footrprint veto
+  muons = "linkedObjects:muons",
+  electrons = "linkedObjects:electrons",
+)
+
+fsrTable = cms.EDProducer("SimpleCandidateFlatTableProducer",
+    src = cms.InputTag("leptonFSRphotons"),
+    cut = cms.string(""), #we should not filter on cross linked collections
+    name = cms.string("FsrPhoton"),
+    doc  = cms.string("Final state radiation photons emitted by muons or electrons"),
+    singleton = cms.bool(False), # the number of entries is variable
+    extension = cms.bool(False), # this is the main table for the muons
+    variables = cms.PSet(P3Vars,
+        relIso03 = Var("userFloat('relIso03')",float,doc="relative isolation in a 0.3 cone without CHS"),
+        dROverEt2 = Var("userFloat('dROverEt2')",float,doc="deltaR to associated muon divided by photon et2"),
+        muonIdx = Var("?hasUserCand('associatedMuon')?userCand('associatedMuon').key():-1",int, doc="index of associated muon"),
+        electronIdx = Var("?hasUserCand('associatedElectron')?userCand('associatedElectron').key():-1",int, doc="index of associated electron")
+        )
+    )
+
+fsrTablesTask =  cms.Task(leptonFSRphotons,fsrTable)

PhysicsTools/NanoAOD/python/muons_cff.py

@@ -87,32 +87,6 @@
         weightFile = "PhysicsTools/NanoAOD/data/mu_BDTG_2016.weights.xml",
     )
 
-from MuonAnalysis.MuonAssociators.muonFSRProducer_cfi import muonFSRProducer
-muonFSRphotons = muonFSRProducer.clone(
-  packedPFCandidates = cms.InputTag("packedPFCandidates"),
-  slimmedElectrons = cms.InputTag("slimmedElectrons"),
-  muons = cms.InputTag("linkedObjects","muons"),
-)
-from MuonAnalysis.MuonAssociators.muonFSRAssociator_cfi import muonFSRAssociator
-muonFSRassociation = muonFSRAssociator.clone(
-    photons = cms.InputTag("muonFSRphotons"),
-    muons = cms.InputTag("linkedObjects","muons"),
-)
-
-fsrTable = cms.EDProducer("SimpleCandidateFlatTableProducer",
-    src = cms.InputTag("muonFSRphotons"),
-    cut = cms.string(""), #we should not filter on cross linked collections
-    name = cms.string("FsrPhoton"),
-    doc  = cms.string("Final state radiation photons emitted by muons"),
-    singleton = cms.bool(False), # the number of entries is variable
-    extension = cms.bool(False), # this is the main table for the muons
-    variables = cms.PSet(P3Vars,
-        relIso03 = Var("userFloat('relIso03')",float,doc="relative isolation in a 0.3 cone without CHS"),
-        dROverEt2 = Var("userFloat('dROverEt2')",float,doc="deltaR to associated muon divided by photon et2"),
-        muonIdx = Var("?hasUserCand('associatedMuon')?userCand('associatedMuon').key():-1",int, doc="index of associated muon")
-        )
-    )
-
 muonTable = cms.EDProducer("SimpleCandidateFlatTableProducer",
     src = cms.InputTag("linkedObjects","muons"),
     cut = cms.string(""), #we should not filter on cross linked collections
@@ -170,7 +144,7 @@
     externalVariables = cms.PSet(
         mvaTTH = ExtVar(cms.InputTag("muonMVATTH"),float, doc="TTH MVA lepton ID score",precision=14),
         mvaLowPt = ExtVar(cms.InputTag("muonMVALowPt"),float, doc="Low pt muon ID score",precision=14),
-        fsrPhotonIdx = ExtVar(cms.InputTag("muonFSRassociation:fsrIndex"),int, doc="Index of the associated FSR photon"),
+        fsrPhotonIdx = ExtVar(cms.InputTag("leptonFSRphotons:muFsrIndex"),int, doc="Index of the lowest-dR/ET2 among associated FSR photons"),
     ),
 )
 
@@ -204,4 +178,4 @@
 
 muonTask = cms.Task(slimmedMuonsUpdated,isoForMu,ptRatioRelForMu,slimmedMuonsWithUserData,finalMuons,finalLooseMuons )
 muonMCTask = cms.Task(muonsMCMatchForTable,muonMCTable)
-muonTablesTask = cms.Task(muonFSRphotons,muonFSRassociation,muonMVATTH,muonMVALowPt,muonTable,fsrTable)
+muonTablesTask = cms.Task(muonMVATTH,muonMVALowPt,muonTable)

PhysicsTools/NanoAOD/python/nano_cff.py

@@ -23,7 +23,7 @@
 from PhysicsTools.NanoAOD.protons_cff import *
 from PhysicsTools.NanoAOD.btagWeightTable_cff import *
 from PhysicsTools.NanoAOD.NanoAODEDMEventContent_cff import *
-
+from PhysicsTools.NanoAOD.fsrPhotons_cff import *
 
 nanoMetadata = cms.EDProducer("UniqueStringProducer",
     strings = cms.PSet(
@@ -66,7 +66,7 @@
      electronTask , lowPtElectronTask, photonTask,
      vertexTask, isoTrackTask, jetLepTask,  # must be after all the leptons
      cms.Task(linkedObjects),
-     jetTablesTask, muonTablesTask, tauTablesTask, boostedTauTablesTask,
+     jetTablesTask, muonTablesTask, fsrTablesTask, tauTablesTask, boostedTauTablesTask,
      electronTablesTask, lowPtElectronTablesTask, photonTablesTask,
      globalTablesTask, vertexTablesTask, metTablesTask, simpleCleanerTable, extraFlagsTableTask,
      isoTrackTablesTask