Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Include FSR for electrons (in addition to muons) #36323

Merged
merged 5 commits into from Jan 11, 2022
Merged
Show file tree
Hide file tree
Changes from 4 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
286 changes: 286 additions & 0 deletions 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& isoConeMax,
const double& isoConeMin) 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, 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, 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& isoConeMax,
const double& isoConeMin) const {
double ptsum = 0;

for (const auto& pfcand : pfcands) {
// Isolation cone
double dRIsoCone = deltaR(photon.eta(), photon.phi(), pfcand.eta(), pfcand.phi());
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

deltaR2 would be computationally cheaper

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Indeed but that would make code less readable, and the same comment applies to the class we are replacing (MuonFSRProducer/MuonFSRAssociator). We can do that if you insist but I have a preference in keeping it more readable.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Uhm... I did not make the same comment in the other place where deltaR is computed exactly because it could have been less readable there, but here dRIsoCone is only used in the comparisons: if you pass (e.g.) isoConeMax2 and isoConeMin2 (and modify accordingly L250) it would be equally readable, and one square root per candidate and per lepton can be avoided

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

OK, done

if (dRIsoCone > isoConeMax || dRIsoCone < isoConeMin)
continue;

// Charged hadrons
if (pfcand.charge() != 0 && abs(pfcand.pdgId()) == 211 && pfcand.pt() > 0.2 && dRIsoCone > drSafe) {
ptsum += pfcand.pt();
// Neutral hadrons + photons
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please indent correctly

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Such indentation is forced by code-format, can't do anything about it.

} else if (pfcand.charge() == 0 && (abs(pfcand.pdgId()) == 22 || abs(pfcand.pdgId()) == 130) && pfcand.pt() > 0.5 &&
dRIsoCone > 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);
4 changes: 3 additions & 1 deletion PhysicsTools/NanoAOD/python/electrons_cff.py
Expand Up @@ -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"),
),
)

Expand Down Expand Up @@ -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


Expand Down
27 changes: 27 additions & 0 deletions 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)
30 changes: 2 additions & 28 deletions PhysicsTools/NanoAOD/python/muons_cff.py
Expand Up @@ -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
Expand Down Expand Up @@ -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"),
),
)

Expand Down Expand Up @@ -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)
4 changes: 2 additions & 2 deletions PhysicsTools/NanoAOD/python/nano_cff.py
Expand Up @@ -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(
Expand Down Expand Up @@ -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
Expand Down