[14_1_X] Adding unbiased superclusters to miniAOD and EGM nano

PhysicsTools/NanoAOD/plugins/SimpleFlatTableProducerPlugins.cc

@@ -6,6 +6,9 @@ typedef SimpleFlatTableProducer<reco::Candidate> SimpleCandidateFlatTableProduce
 #include "DataFormats/TrackReco/interface/Track.h"
 typedef SimpleFlatTableProducer<reco::Track> SimpleTrackFlatTableProducer;
 
+#include "DataFormats/EgammaReco/interface/SuperCluster.h"
+typedef SimpleFlatTableProducer<reco::SuperCluster> SimpleSuperclusterFlatTableProducer;
+
 #include "DataFormats/JetReco/interface/PFJet.h"
 typedef SimpleFlatTableProducer<reco::PFJet> SimplePFJetFlatTableProducer;
 
@@ -48,6 +51,7 @@ typedef EventSingletonSimpleFlatTableProducer<reco::BeamSpot> SimpleBeamspotFlat
 #include "FWCore/Framework/interface/MakerMacros.h"
 DEFINE_FWK_MODULE(SimpleCandidateFlatTableProducer);
 DEFINE_FWK_MODULE(SimpleTrackFlatTableProducer);
+DEFINE_FWK_MODULE(SimpleSuperclusterFlatTableProducer);
 DEFINE_FWK_MODULE(SimplePFJetFlatTableProducer);
 DEFINE_FWK_MODULE(SimpleVertexFlatTableProducer);
 DEFINE_FWK_MODULE(SimpleSecondaryVertexFlatTableProducer);

PhysicsTools/NanoAOD/python/egamma_custom_cff.py

@@ -8,6 +8,7 @@
 from PhysicsTools.NanoAOD.nanoDQM_cfi import nanoDQM
 from PhysicsTools.NanoAOD.nanoDQM_cff import _Photon_extra_plots, _Electron_extra_plots
 from PhysicsTools.NanoAOD.triggerObjects_cff import triggerObjectTable, mksel
+from RecoEgamma.EgammaIsolationAlgos.superclusValueMapProducer_cfi import superclusValueMaps
 
 customElectronFilterBits = cms.PSet(
     doc = cms.string("PixelMatched e/gamma"), # this may also select photons!
@@ -81,6 +82,27 @@
     )
 )
 
+superclusterTable = cms.EDProducer("SimpleSuperclusterFlatTableProducer",
+  src = cms.InputTag("reducedEgamma","reducedSuperClusters"),
+  name = cms.string("Supercluster"),
+  doc = cms.string("Supercluster collection"),
+  variables = cms.PSet(
+    energy = Var("energy()",float,doc="supercluster energy",precision=10),
+    eta = Var("eta()",float,doc="supercluster eta",precision=10),
+    phi = Var("phi()",float,doc="supercluster phi",precision=10),
+    rawEnergy = Var("rawEnergy()",float,doc="sum of basic clusters energy",precision=10),
+    preshowerEnergy = Var("preshowerEnergy()",float,doc="sum of energy in preshower",precision=10),
+    etaWidth = Var("etaWidth()",float,doc="supercluster eta width",precision=10),
+    phiWidth = Var("etaWidth()",float,doc="supercluster phi width",precision=10),
+    seedClusEnergy = Var("seed().energy()",float,doc="seed cluster energy",precision=10),
+    seedClusterEta = Var("seed().eta()",float,doc="seed cluster eta",precision=10),
+    seedClusterPhi = Var("seed().phi()",float,doc="seed cluster phi",precision=10),
+  ),
+  externalVariables = cms.PSet(
+    trkIso = ExtVar("superclusValueMaps:superclusTkIso",float,doc="supercluster track iso within 0.06 < dR < 0.4 & |dEta| > 0.03",precision=10),
+  )
+)
+
 def addExtraEGammaVarsCustomize(process):
     #photon
     process.photonTable.variables.setValue(_phoVarsExtra.parameters_())
@@ -92,5 +114,13 @@ def addExtraEGammaVarsCustomize(process):
     process.triggerObjectTable.selections.Electron = customElectronFilterBits
     if process.nanoDQM:
       process.nanoDQM.vplots.Electron.plots = _Electron_extra_plots
+
+    #superCluster
+    process.superclusValueMaps = superclusValueMaps
+    process.superclusterTable = superclusterTable
+
+    process.superclusterTask = cms.Task(process.superclusValueMaps)
+    process.superclusterTask.add(process.superclusterTable)
+    process.nanoTableTaskCommon.add(process.superclusterTask)
 
     return process

PhysicsTools/NanoAOD/python/electrons_cff.py

@@ -431,7 +431,7 @@ def _get_bitmapVIDForEle_docstring(modules,WorkingPoints):
         electronTable.variables,
         pt = Var("pt*userFloat('ecalTrkEnergyPostCorrNew')/userFloat('ecalTrkEnergyPreCorrNew')", float, precision=-1, doc="p_{T}"),
         energyErr = Var("userFloat('ecalTrkEnergyErrPostCorrNew')", float, precision=6, doc="energy error of the cluster-track combination"),
-        eCorr = Var("userFloat('ecalTrkEnergyPostCorrNew')/userFloat('ecalTrkEnergyPreCorrNew')", float, doc="ratio of the calibrated energy/miniaod energy"),
+        ptPreCorr = Var("pt", float, doc="pt of the electron before energy corrections"),
         scEtOverPt = Var("(superCluster().energy()/(pt*userFloat('ecalTrkEnergyPostCorrNew')/userFloat('ecalTrkEnergyPreCorrNew')*cosh(superCluster().eta())))-1",float,doc="(supercluster transverse energy)/pt-1",precision=8),
         dEscaleUp=Var("userFloat('ecalTrkEnergyPostCorrNew')-userFloat('energyScaleUpNew')", float,  doc="ecal energy scale shifted 1 sigma up(adding gain/stat/syst in quadrature)", precision=8),
         dEscaleDown=Var("userFloat('ecalTrkEnergyPostCorrNew')-userFloat('energyScaleDownNew')", float,  doc="ecal energy scale shifted 1 sigma down (adding gain/stat/syst in quadrature)", precision=8),

PhysicsTools/NanoAOD/python/nanoDQM_cff.py

@@ -123,7 +123,7 @@
     Plot1D('dEscaleDown', 'dEscaleDown', 100, -0.01, 0.01, '#Delta E scaleDown'),
     Plot1D('dEsigmaUp', 'dEsigmaUp', 100, -0.1, 0.1, '#Delta E sigmaUp'),
     Plot1D('dEsigmaDown', 'dEsigmaDown', 100, -0.1, 0.1, '#Delta E sigmaDown'),
-    Plot1D('eCorr', 'eCorr', 20, 0.8, 1.2, 'ratio of the calibrated energy/miniaod energy'),
+    Plot1D('ptPreCorr', 'ptPreCorr', 100, 0., 500., 'Pt before scale & smearing energy corrections'),
 ])
 run2_egamma.toModify(
      nanoDQM.vplots.Electron,
@@ -146,7 +146,7 @@ def _match(name):
     Plot1D('dEscaleDown', 'dEscaleDown', 100, -0.01, 0.01, '#Delta E scaleDown'),
     Plot1D('dEsigmaUp', 'dEsigmaUp', 100, -0.1, 0.1, '#Delta E sigmaUp'),
     Plot1D('dEsigmaDown', 'dEsigmaDown', 100, -0.1, 0.1, '#Delta E sigmaDown'),
-    Plot1D('eCorr', 'eCorr', 20, 0.8, 1.2, 'ratio of the calibrated energy/miniaod energy'),
+    Plot1D('ptPreCorr', 'ptPreCorr', 100, 0., 500., 'Pt before scale & smearing energy corrections'),
 ])
 run2_egamma.toModify(
      nanoDQM.vplots.Photon,

PhysicsTools/NanoAOD/python/photons_cff.py

@@ -283,7 +283,7 @@ def make_bitmapVID_docstring(id_modules_working_points_pset):
     photonTable.variables,
     pt = Var("pt*userFloat('ecalEnergyPostCorrNew')/userFloat('ecalEnergyPreCorrNew')", float, precision=-1, doc="p_{T}"),
     energyErr = Var("userFloat('ecalEnergyErrPostCorrNew')",float,doc="energy error of the cluster from regression",precision=6),
-    eCorr = Var("userFloat('ecalEnergyPostCorrNew')/userFloat('ecalEnergyPreCorrNew')",float,doc="ratio of the calibrated energy/miniaod energy"),
+    ptPreCorr = Var("pt",float,doc="pt of the photon before energy corrections"),
     cutBased = Var(
             "userInt('cutBasedID_Fall17V2_loose')+userInt('cutBasedID_Fall17V2_medium')+userInt('cutBasedID_Fall17V2_tight')",
             "uint8",

PhysicsTools/PatAlgos/python/slimming/slimming_cff.py

@@ -66,6 +66,7 @@
     slimmedLambdaVertices,
     slimmedMETs,
     metFilterPathsTask,
+    superClusterMerger,
     reducedEgamma,
     slimmedHcalRecHits,
     bunchSpacingProducer,

RecoEgamma/EgammaIsolationAlgos/interface/EleTkIsolFromCands.h

@@ -92,6 +92,12 @@ class EleTkIsolFromCands {
 
   Output operator()(const reco::TrackBase& electronTrack);
 
+protected:
+  using TrackTable = edm::soa::Table<edm::soa::col::Pt, edm::soa::col::Eta, edm::soa::col::Phi, edm::soa::col::Vz>;
+  TrackTable const& getPreselectedTracks(bool isBarrel);
+
+  Configuration const& cfg_;
+
 private:
   // For each electron, we want to try out which tracks are in a cone around
   // it. However, this will get expensive if there are many electrons and
@@ -103,8 +109,6 @@ class EleTkIsolFromCands {
   // the electron. Note that this has to be done twice, because the required
   // preselection is different for barrel and endcap electrons.
 
-  using TrackTable = edm::soa::Table<edm::soa::col::Pt, edm::soa::col::Eta, edm::soa::col::Phi, edm::soa::col::Vz>;
-
   static bool passPIDVeto(const int pdgId, const EleTkIsolFromCands::PIDVeto pidVeto);
 
   static TrackTable preselectTracks(reco::TrackCollection const& tracks, TrkCuts const& cuts);
@@ -118,10 +122,6 @@ class EleTkIsolFromCands {
   static bool passQual(const reco::TrackBase& trk, const std::vector<reco::TrackBase::TrackQuality>& quals);
   static bool passAlgo(const reco::TrackBase& trk, const std::vector<reco::TrackBase::TrackAlgorithm>& algosToRej);
 
-  TrackTable const& getPreselectedTracks(bool isBarrel);
-
-  Configuration const& cfg_;
-
   // All of these member variables are related to the caching of preselected tracks
   reco::TrackCollection const* tracks_ = nullptr;
   pat::PackedCandidateCollection const* cands_ = nullptr;

RecoEgamma/EgammaIsolationAlgos/interface/SuperclusTkIsolFromCands.h

@@ -0,0 +1,27 @@
+#ifndef RecoEgamma_EgammaIsolationAlgos_SuperclusTkIsolFromCands_h
+#define RecoEgamma_EgammaIsolationAlgos_SuperclusTkIsolFromCands_h
+
+#include "CommonTools/Utils/interface/KinematicColumns.h"
+#include "DataFormats/TrackReco/interface/TrackBase.h"
+#include "DataFormats/TrackReco/interface/TrackFwd.h"
+#include "DataFormats/PatCandidates/interface/PackedCandidate.h"
+#include "DataFormats/Math/interface/Point3D.h"
+#include "DataFormats/EgammaReco/interface/SuperCluster.h"
+#include "FWCore/SOA/interface/Table.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
+#include "RecoEgamma/EgammaIsolationAlgos/interface/EleTkIsolFromCands.h"
+
+class SuperclusTkIsolFromCands : public EleTkIsolFromCands {
+public:
+  explicit SuperclusTkIsolFromCands(Configuration const& cfg, reco::TrackCollection const& tracks)
+      : EleTkIsolFromCands(cfg, tracks) {}
+  explicit SuperclusTkIsolFromCands(Configuration const& cfg,
+                                    pat::PackedCandidateCollection const& cands,
+                                    PIDVeto pidVeto = PIDVeto::NONE)
+      : EleTkIsolFromCands(cfg, cands, pidVeto) {}
+
+  Output operator()(const reco::SuperCluster& sc, const math::XYZPoint& vtx);
+};
+
+#endif

RecoEgamma/EgammaIsolationAlgos/plugins/SuperclusValueMapProducer.cc

@@ -0,0 +1,129 @@
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+#include "FWCore/Framework/interface/stream/EDProducer.h"
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/ESHandle.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+
+#include "DataFormats/Common/interface/ValueMap.h"
+#include "DataFormats/Common/interface/View.h"
+#include "DataFormats/Common/interface/Handle.h"
+#include "DataFormats/Math/interface/Point3D.h"
+#include "DataFormats/PatCandidates/interface/PackedCandidate.h"
+
+#include "RecoEgamma/EgammaIsolationAlgos/interface/SuperclusTkIsolFromCands.h"
+
+// copy-pasted from ElectronHEEPIDValueMapProducer
+
+class SuperclusValueMapProducer : public edm::stream::EDProducer<> {
+public:
+  explicit SuperclusValueMapProducer(const edm::ParameterSet&);
+  ~SuperclusValueMapProducer() = default;
+
+private:
+  void produce(edm::Event&, const edm::EventSetup&) override;
+
+  std::vector<edm::EDGetTokenT<pat::PackedCandidateCollection>> setTokens(const std::vector<edm::InputTag>& tags);
+
+  template <typename T>
+  static void writeValueMap(edm::Event& iEvent,
+                            const edm::Handle<edm::View<reco::SuperCluster>>& handle,
+                            const std::vector<T>& values,
+                            const std::string& label);
+
+  std::vector<SuperclusTkIsolFromCands::PIDVeto> candVetos_;
+
+  const std::vector<edm::InputTag> candTags_;
+  const std::vector<edm::EDGetTokenT<pat::PackedCandidateCollection>> candTokens_;
+  const edm::EDGetTokenT<edm::View<reco::SuperCluster>> scToken_;
+  const edm::EDGetTokenT<edm::View<reco::Vertex>> pvToken_;
+  const edm::EDGetTokenT<reco::BeamSpot> bsToken_;
+
+  const SuperclusTkIsolFromCands::Configuration trkIsoCalcCfg_;
+
+  const std::string superclusTkIsoLabel_ = "superclusTkIso";
+};
+
+SuperclusValueMapProducer::SuperclusValueMapProducer(const edm::ParameterSet& iConfig)
+    : candTags_(iConfig.getParameter<std::vector<edm::InputTag>>("cands")),
+      candTokens_(setTokens(candTags_)),
+      scToken_(consumes<edm::View<reco::SuperCluster>>(iConfig.getParameter<edm::InputTag>("srcSc"))),
+      pvToken_(consumes<edm::View<reco::Vertex>>(iConfig.getParameter<edm::InputTag>("srcPv"))),
+      bsToken_(consumes<reco::BeamSpot>(iConfig.getParameter<edm::InputTag>("srcBs"))),
+      trkIsoCalcCfg_(iConfig.getParameter<edm::ParameterSet>("trkIsoConfig")) {
+  auto fillVetos = [](const auto& in, auto& out) {
+    std::transform(in.begin(), in.end(), std::back_inserter(out), SuperclusTkIsolFromCands::pidVetoFromStr);
+  };
+
+  fillVetos(iConfig.getParameter<std::vector<std::string>>("candVetos"), candVetos_);
+
+  if (candVetos_.size() != candTags_.size())
+    throw cms::Exception("ConfigError") << "Error candVetos should be the same size as cands" << std::endl;
+
+  produces<edm::ValueMap<float>>(superclusTkIsoLabel_);
+}
+
+void SuperclusValueMapProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
+  edm::Handle<edm::View<reco::SuperCluster>> scHandle;
+  iEvent.getByToken(scToken_, scHandle);
+
+  edm::Handle<edm::View<reco::Vertex>> pvHandle;
+  iEvent.getByToken(pvToken_, pvHandle);
+
+  edm::Handle<reco::BeamSpot> bsHandle;
+  iEvent.getByToken(bsToken_, bsHandle);
+
+  math::XYZPoint pos;
+
+  if (pvHandle.isValid() && !pvHandle->empty())
+    pos = pvHandle->front().position();  // first try PV
+  else
+    pos = (*bsHandle).position();  // fall back to BS
+
+  std::vector<edm::Handle<pat::PackedCandidateCollection>> candHandles(candTokens_.size());
+  std::vector<std::unique_ptr<SuperclusTkIsolFromCands>> tkIsoCalc;
+
+  for (unsigned idx = 0; idx < candTokens_.size(); idx++) {
+    iEvent.getByToken(candTokens_.at(idx), candHandles.at(idx));
+    tkIsoCalc.push_back(
+        std::make_unique<SuperclusTkIsolFromCands>(trkIsoCalcCfg_, *(candHandles.at(idx)), candVetos_.at(idx)));
+  }
+
+  std::vector<float> vecTkIso;
+  vecTkIso.reserve(scHandle->size());
+
+  for (const auto& sc : *scHandle) {
+    float tkIso = 0.;
+
+    for (auto& calc : tkIsoCalc)
+      tkIso += (*calc)(sc, pos).ptSum;
+
+    vecTkIso.push_back(tkIso);
+  }
+
+  writeValueMap(iEvent, scHandle, vecTkIso, superclusTkIsoLabel_);
+}
+
+std::vector<edm::EDGetTokenT<pat::PackedCandidateCollection>> SuperclusValueMapProducer::setTokens(
+    const std::vector<edm::InputTag>& tags) {
+  std::vector<edm::EDGetTokenT<pat::PackedCandidateCollection>> out;
+
+  for (const auto& tag : tags)
+    out.push_back(consumes<pat::PackedCandidateCollection>(tag));
+
+  return out;
+}
+
+template <typename T>
+void SuperclusValueMapProducer::writeValueMap(edm::Event& iEvent,
+                                              const edm::Handle<edm::View<reco::SuperCluster>>& handle,
+                                              const std::vector<T>& values,
+                                              const std::string& label) {
+  std::unique_ptr<edm::ValueMap<T>> valMap(new edm::ValueMap<T>());
+  typename edm::ValueMap<T>::Filler filler(*valMap);
+  filler.insert(handle, values.begin(), values.end());
+  filler.fill();
+  iEvent.put(std::move(valMap), label);
+}
+
+DEFINE_FWK_MODULE(SuperclusValueMapProducer);

RecoEgamma/EgammaIsolationAlgos/python/superclusValueMapProducer_cfi.py

@@ -0,0 +1,33 @@
+import FWCore.ParameterSet.Config as cms
+
+from RecoEgamma.EgammaIsolationAlgos.electronTrackIsolations_cfi import trkIsol04CfgV2
+
+# supposed to be calculated from AOD (see https://github.com/cms-sw/cmssw/pull/42007)
+# but since we didn't create the new PAT version of SC
+# trk iso values are calculated on-the-fly from miniAOD
+# parameters are inherited from HEEP trk iso
+
+# instead, we use fairly loose inner cone & strip veto from EgammaHLTTrackIsolation
+# to avoid strong bias due to the tight trk isolation
+# e.g. see hltEgammaHollowTrackIsoL1Seeded filter
+
+scTrkIso04 = cms.PSet(
+  barrelCuts = trkIsol04CfgV2.barrelCuts.clone(
+    minDR = 0.06,
+    minDEta = 0.03
+  ),
+  endcapCuts = trkIsol04CfgV2.endcapCuts.clone(
+    minDR = 0.06,
+    minDEta = 0.03
+  )
+)
+
+superclusValueMaps = cms.EDProducer("SuperclusValueMapProducer",
+  srcBs = cms.InputTag("offlineBeamSpot"),
+  srcPv = cms.InputTag("offlineSlimmedPrimaryVertices"),
+  srcSc = cms.InputTag("reducedEgamma:reducedSuperClusters"),
+  cands = cms.VInputTag("packedPFCandidates",
+                        "lostTracks"), # do not count electron tracks in the trk iso
+  candVetos = cms.vstring("ELES","NONE"),
+  trkIsoConfig = scTrkIso04,
+)

RecoEgamma/EgammaIsolationAlgos/src/SuperclusTkIsolFromCands.cc

@@ -0,0 +1,33 @@
+#include "RecoEgamma/EgammaIsolationAlgos/interface/SuperclusTkIsolFromCands.h"
+#include "DataFormats/TrackReco/interface/Track.h"
+#include "DataFormats/Math/interface/deltaR.h"
+
+SuperclusTkIsolFromCands::Output SuperclusTkIsolFromCands::operator()(const reco::SuperCluster& sc,
+                                                                      const math::XYZPoint& vtx) {
+  using namespace edm::soa::col;
+
+  float ptSum = 0.;
+  int nrTrks = 0;
+
+  const float scEta = sc.eta();
+  const float scPhi = sc.phi();
+  const float vtxVz = vtx.z();
+
+  const bool isBarrelSC = std::abs(scEta) < 1.5;
+
+  auto const& preselectedTracks = getPreselectedTracks(isBarrelSC);
+  auto const& cuts = isBarrelSC ? cfg_.barrelCuts : cfg_.endcapCuts;
+
+  for (auto const& trk : preselectedTracks) {
+    const float dR2 = reco::deltaR2(scEta, scPhi, trk.get<Eta>(), trk.get<Phi>());
+    const float dEta = trk.get<Eta>() - scEta;
+    const float dZ = vtxVz - trk.get<Vz>();
+
+    if (dR2 >= cuts.minDR2 && dR2 <= cuts.maxDR2 && std::abs(dEta) >= cuts.minDEta && std::abs(dZ) < cuts.maxDZ) {
+      ptSum += trk.get<Pt>();
+      nrTrks++;
+    }
+  }
+
+  return {nrTrks, ptSum};
+}