Merge pull request cms-sw#91 from cerminar/l2EgFix_v1

Fixes and improvements for Layer-2 e/g objects and quality bits in the endcap

L1Trigger/L1THGCal/interface/backend/HGCalTriggerClusterIdentificationBase.h

@@ -10,7 +10,8 @@ class HGCalTriggerClusterIdentificationBase {
   virtual ~HGCalTriggerClusterIdentificationBase(){};
   virtual void initialize(const edm::ParameterSet& conf) = 0;
   virtual float value(const l1t::HGCalMulticluster& cluster) const = 0;
-  virtual bool decision(const l1t::HGCalMulticluster& cluster) const = 0;
+  virtual bool decision(const l1t::HGCalMulticluster& cluster, unsigned wp = 0) const = 0;
+  virtual const std::vector<std::string>& working_points() const = 0;
 };
 
 #include "FWCore/PluginManager/interface/PluginFactory.h"

L1Trigger/L1THGCal/plugins/backend/HGCalTriggerClusterIdentificationBDT.cc

@@ -33,12 +33,26 @@ class HGCalTriggerClusterIdentificationBDT : public HGCalTriggerClusterIdentific
     float eta_max_ = 3.;
   };
 
+  class WorkingPoint {
+  public:
+    WorkingPoint(const std::string& name, double wp) : name_(name), wp_(wp) {}
+    ~WorkingPoint() = default;
+
+    const std::string& name() const { return name_; }
+    double working_point() const { return wp_; }
+
+  private:
+    std::string name_;
+    double wp_;
+  };
+
 public:
   HGCalTriggerClusterIdentificationBDT();
   ~HGCalTriggerClusterIdentificationBDT() override{};
   void initialize(const edm::ParameterSet& conf) final;
   float value(const l1t::HGCalMulticluster& cluster) const final;
-  bool decision(const l1t::HGCalMulticluster& cluster) const final;
+  bool decision(const l1t::HGCalMulticluster& cluster, unsigned wp = 0) const final;
+  const std::vector<std::string>& working_points() const final { return working_points_names_; }
 
 private:
   enum class ClusterVariable {
@@ -57,9 +71,10 @@ class HGCalTriggerClusterIdentificationBDT : public HGCalTriggerClusterIdentific
   };
   std::vector<Category> categories_;
   std::vector<std::unique_ptr<TMVAEvaluator>> bdts_;
-  std::vector<double> working_points_;
+  std::vector<std::vector<WorkingPoint>> working_points_;
   std::vector<std::string> input_variables_;
   std::vector<ClusterVariable> input_variables_id_;
+  std::vector<std::string> working_points_names_;
 
   float clusterVariable(ClusterVariable, const l1t::HGCalMulticluster&) const;
   int category(float pt, float eta) const;
@@ -81,17 +96,32 @@ void HGCalTriggerClusterIdentificationBDT::initialize(const edm::ParameterSet& c
   std::vector<double> categories_etamax = conf.getParameter<std::vector<double>>("CategoriesEtaMax");
   std::vector<double> categories_ptmin = conf.getParameter<std::vector<double>>("CategoriesPtMin");
   std::vector<double> categories_ptmax = conf.getParameter<std::vector<double>>("CategoriesPtMax");
-  working_points_ = conf.getParameter<std::vector<double>>("WorkingPoints");
 
   if (bdt_files.size() != categories_etamin.size() || categories_etamin.size() != categories_etamax.size() ||
-      categories_etamax.size() != categories_ptmin.size() || categories_ptmin.size() != categories_ptmax.size() ||
-      categories_ptmax.size() != working_points_.size()) {
+      categories_etamax.size() != categories_ptmin.size() || categories_ptmin.size() != categories_ptmax.size()) {
     throw cms::Exception("HGCalTriggerClusterIdentificationBDT|BadInitialization")
         << "Inconsistent numbers of categories, BDT weight files and working points";
   }
-  categories_.reserve(working_points_.size());
-  bdts_.reserve(working_points_.size());
-  for (unsigned cat = 0; cat < categories_etamin.size(); cat++) {
+  size_t categories_size = categories_etamin.size();
+
+  const auto wps_conf = conf.getParameter<std::vector<edm::ParameterSet>>("WorkingPoints");
+  working_points_.resize(categories_size);
+  for (const auto& wp_conf : wps_conf) {
+    std::string wp_name = wp_conf.getParameter<std::string>("Name");
+    std::vector<double> wps = wp_conf.getParameter<std::vector<double>>("WorkingPoint");
+    working_points_names_.emplace_back(wp_name);
+    if (wps.size() != categories_size) {
+      throw cms::Exception("HGCalTriggerClusterIdentificationBDT|BadInitialization")
+          << "Inconsistent number of categories in working point '" << wp_name << "'";
+    }
+    for (size_t cat = 0; cat < categories_size; cat++) {
+      working_points_[cat].emplace_back(wp_name, wps[cat]);
+    }
+  }
+
+  categories_.reserve(categories_size);
+  bdts_.reserve(categories_size);
+  for (size_t cat = 0; cat < categories_size; cat++) {
     categories_.emplace_back(
         categories_ptmin[cat], categories_ptmax[cat], categories_etamin[cat], categories_etamax[cat]);
   }
@@ -148,12 +178,12 @@ float HGCalTriggerClusterIdentificationBDT::value(const l1t::HGCalMulticluster&
   return (cat != -1 ? bdts_.at(cat)->evaluate(inputs) : -999.);
 }
 
-bool HGCalTriggerClusterIdentificationBDT::decision(const l1t::HGCalMulticluster& cluster) const {
+bool HGCalTriggerClusterIdentificationBDT::decision(const l1t::HGCalMulticluster& cluster, unsigned wp) const {
   float bdt_output = value(cluster);
   float pt = cluster.pt();
   float eta = cluster.eta();
   int cat = category(pt, eta);
-  return (cat != -1 ? bdt_output > working_points_.at(cat) : true);
+  return (cat != -1 ? bdt_output > working_points_.at(cat).at(wp).working_point() : true);
 }
 
 int HGCalTriggerClusterIdentificationBDT::category(float pt, float eta) const {

L1Trigger/L1THGCal/python/egammaIdentification.py

@@ -205,23 +205,44 @@ def __init__(self, eta_min, eta_max, pt_min, pt_max):
         CategoriesPtMin=cms.vdouble([cat.pt_min for cat in categories]),
         CategoriesPtMax=cms.vdouble([cat.pt_max for cat in categories]),
         Weights=cms.vstring(bdt_weights_histomax['v10_3151']),
-        WorkingPoints=cms.vdouble(
-                [wps[eff] for wps,eff in zip(working_points_histomax['v10_3151'],tight_wp)]
-                )
+        WorkingPoints=cms.VPSet([
+            cms.PSet(
+                Name=cms.string('tight'),
+                WorkingPoint=cms.vdouble([wps[eff] for wps,eff in zip(working_points_histomax['v10_3151'],tight_wp)])
+            ),
+            cms.PSet(
+                Name=cms.string('loose'),
+                WorkingPoint=cms.vdouble([wps[eff] for wps,eff in zip(working_points_histomax['v10_3151'],loose_wp)])
+            ),
+            ])
         )
 
 phase2_hgcalV10.toModify(egamma_identification_histomax,
         Inputs=cms.vstring(input_features_histomax['v10_3151']),
         Weights=cms.vstring(bdt_weights_histomax['v10_3151']),
-        WorkingPoints=cms.vdouble(
-                [wps[eff] for wps,eff in zip(working_points_histomax['v10_3151'],tight_wp)]
-                )
+        WorkingPoints=cms.VPSet([
+            cms.PSet(
+                Name=cms.string('tight'),
+                WorkingPoint=cms.vdouble([wps[eff] for wps,eff in zip(working_points_histomax['v10_3151'],tight_wp)])
+            ),
+            cms.PSet(
+                Name=cms.string('loose'),
+                WorkingPoint=cms.vdouble([wps[eff] for wps,eff in zip(working_points_histomax['v10_3151'],loose_wp)])
+            ),
+            ])
         )
 
 phase2_hgcalV11.toModify(egamma_identification_histomax,
         Inputs=cms.vstring(input_features_histomax['v10_3151']),
         Weights=cms.vstring(bdt_weights_histomax['v10_3151']),
-        WorkingPoints=cms.vdouble(
-                [wps[eff] for wps,eff in zip(working_points_histomax['v10_3151'],tight_wp)]
-                )
+        WorkingPoints=cms.VPSet([
+            cms.PSet(
+                Name=cms.string('tight'),
+                WorkingPoint=cms.vdouble([wps[eff] for wps,eff in zip(working_points_histomax['v10_3151'],tight_wp)])
+            ),
+            cms.PSet(
+                Name=cms.string('loose'),
+                WorkingPoint=cms.vdouble([wps[eff] for wps,eff in zip(working_points_histomax['v10_3151'],loose_wp)])
+            ),
+            ])
         )

L1Trigger/L1THGCal/src/backend/HGCalHistoClusteringImpl.cc

@@ -141,7 +141,11 @@ void HGCalHistoClusteringImpl::finalizeClusters(std::vector<l1t::HGCalMulticlust
       //compute shower shapes
       shape_.fillShapes(multicluster, triggerGeometry);
       // fill quality flag
-      multicluster.setHwQual(id_->decision(multicluster));
+      unsigned hwQual = 0;
+      for (unsigned wp = 0; wp < id_->working_points().size(); wp++) {
+        hwQual |= (id_->decision(multicluster, wp) << wp);
+      }
+      multicluster.setHwQual(hwQual);
       // fill H/E
       multicluster.saveHOverE();
 

L1Trigger/Phase2L1ParticleFlow/plugins/L1TCtL2EgProducer.cc

@@ -49,7 +49,7 @@ class L1TCtL2EgProducer : public edm::global::EDProducer<> {
 
     BXVector<edm::Ref<BXVector<l1t::EGamma>>> oldRefs;
     std::map<edm::Ref<BXVector<l1t::EGamma>>, edm::Ref<BXVector<l1t::EGamma>>> old2newRefMap;
-    std::vector<std::pair<const edm::Ref<l1t::EGammaBxCollection> &, edm::Ptr<L1TTTrackType>>> origRefAndPtr;
+    std::vector<std::pair<edm::Ref<l1t::EGammaBxCollection>, edm::Ptr<L1TTTrackType>>> origRefAndPtr;
   };
 
   void convertToEmu(const l1t::TkElectron &tkele, RefRemapper &refRemapper, l1ct::OutputBoard &boarOut) const;
@@ -362,34 +362,36 @@ void L1TCtL2EgProducer::convertToEmu(const l1t::TkEm &tkem,
 
 l1t::TkEm L1TCtL2EgProducer::convertFromEmu(const l1ct::EGIsoObjEmu &egiso, const RefRemapper &refRemapper) const {
   // std::cout << "[convertFromEmu] TkEm pt: " << egiso.hwPt << " eta: " << egiso.hwEta << " phi: " << egiso.hwPhi << " staidx: " << egiso.sta_idx << std::endl;
-
-  reco::Candidate::PolarLorentzVector mom(egiso.floatPt(), egiso.floatEta(), egiso.floatPhi(), 0.);
+  // NOTE: the TkEM object is created with the accuracy as in GT object (not the Correlator internal one)!
+  const auto gteg = egiso.toGT();
+  reco::Candidate::PolarLorentzVector mom(
+      l1gt::Scales::floatPt(gteg.v3.pt), l1gt::Scales::floatEta(gteg.v3.eta), l1gt::Scales::floatPhi(gteg.v3.phi), 0.);
   l1t::TkEm tkem(reco::Candidate::LorentzVector(mom),
                  refRemapper.origRefAndPtr[egiso.sta_idx].first,
                  egiso.floatRelIso(l1ct::EGIsoObjEmu::IsoType::TkIso),
                  egiso.floatRelIso(l1ct::EGIsoObjEmu::IsoType::TkIsoPV));
-  // FIXME: need to define a global quality (barrel+endcap) or add a bit to distibguish them?
-  tkem.setHwQual(egiso.hwQual);
+  tkem.setHwQual(gteg.quality);
   tkem.setPFIsol(egiso.floatRelIso(l1ct::EGIsoObjEmu::IsoType::PfIso));
   tkem.setPFIsolPV(egiso.floatRelIso(l1ct::EGIsoObjEmu::IsoType::PfIsoPV));
-  tkem.setEgBinaryWord(egiso.toGT().pack());
+  tkem.setEgBinaryWord(gteg.pack());
   return tkem;
 }
 
 l1t::TkElectron L1TCtL2EgProducer::convertFromEmu(const l1ct::EGIsoEleObjEmu &egele,
                                                   const RefRemapper &refRemapper) const {
   // std::cout << "[convertFromEmu] TkEle pt: " << egele.hwPt << " eta: " << egele.hwEta << " phi: " << egele.hwPhi << " staidx: " << egele.sta_idx << std::endl;
-
-  reco::Candidate::PolarLorentzVector mom(egele.floatPt(), egele.hwEta, egele.hwPhi, 0.);
+  // NOTE: the TkElectron object is created with the accuracy as in GT object (not the Correlator internal one)!
+  const auto gteg = egele.toGT();
+  reco::Candidate::PolarLorentzVector mom(
+      l1gt::Scales::floatPt(gteg.v3.pt), l1gt::Scales::floatEta(gteg.v3.eta), l1gt::Scales::floatPhi(gteg.v3.phi), 0.);
 
   l1t::TkElectron tkele(reco::Candidate::LorentzVector(mom),
                         refRemapper.origRefAndPtr[egele.sta_idx].first,
                         refRemapper.origRefAndPtr[egele.sta_idx].second,
                         egele.floatRelIso(l1ct::EGIsoEleObjEmu::IsoType::TkIso));
-  // FIXME: need to define a global quality (barrel+endcap)?
-  tkele.setHwQual(egele.hwQual);
+  tkele.setHwQual(gteg.quality);
   tkele.setPFIsol(egele.floatRelIso(l1ct::EGIsoEleObjEmu::IsoType::PfIso));
-  tkele.setEgBinaryWord(egele.toGT().pack());
+  tkele.setEgBinaryWord(gteg.pack());
   return tkele;
 }
 

L1Trigger/Phase2L1ParticleFlow/python/l1TkEgAlgoEmulator_cfi.py

@@ -10,6 +10,7 @@
     writeBeforeBremRecovery=cms.bool(True),
     filterHwQuality=cms.bool(False),
     caloHwQual=cms.int32(4),
+    doEndcapHwQual=cms.bool(False),
     dEtaMaxBrem=cms.double(0.02),
     dPhiMaxBrem=cms.double(0.1),
     absEtaBoundaries=cms.vdouble(0.0, 0.9, 1.5),

L1Trigger/Phase2L1ParticleFlow/python/l1ctLayer1_cff.py

@@ -255,6 +255,7 @@
         nEMCALO_EGIN = 10, 
         nEM_EGOUT = 5,
         doBremRecovery=True,
+        doEndcapHwQual=True,
         writeBeforeBremRecovery=False,
         writeEGSta=True),
     tkEgSorterParameters=tkEgSorterParameters.clone(
@@ -340,6 +341,7 @@
         nEMCALO_EGIN = 10, 
         nEM_EGOUT = 5,
         doBremRecovery=True,
+        doEndcapHwQual=True,
         writeBeforeBremRecovery=False,
         writeEGSta=True),
     tkEgSorterParameters=tkEgSorterParameters.clone(

L1Trigger/Phase2L1ParticleFlow/src/newfirmware/egamma/pftkegalgo_ref.cpp

@@ -5,6 +5,7 @@
 #include <algorithm>
 #include <memory>
 #include <iostream>
+#include <bitset>
 
 #include "L1Trigger/Phase2L1ParticleFlow/src/dbgPrintf.h"
 
@@ -22,6 +23,7 @@ l1ct::PFTkEGAlgoEmuConfig::PFTkEGAlgoEmuConfig(const edm::ParameterSet &pset)
       doBremRecovery(pset.getParameter<bool>("doBremRecovery")),
       writeBeforeBremRecovery(pset.getParameter<bool>("writeBeforeBremRecovery")),
       caloHwQual(pset.getParameter<int>("caloHwQual")),
+      doEndcapHwQual(pset.getParameter<bool>("doEndcapHwQual")),
       emClusterPtMin(pset.getParameter<double>("caloEtMin")),
       dEtaMaxBrem(pset.getParameter<double>("dEtaMaxBrem")),
       dPhiMaxBrem(pset.getParameter<double>("dPhiMaxBrem")),
@@ -166,7 +168,8 @@ void PFTkEGAlgoEmulator::run(const PFInputRegion &in, OutputRegion &out) const {
       if (calo.hwPt > 0)
         dbgCout() << "[REF] IN calo[" << ic << "] pt: " << calo.hwPt << " eta: " << calo.hwEta
                   << " (glb eta: " << in.region.floatGlbEta(calo.hwEta) << ") phi: " << calo.hwPhi
-                  << "(glb phi: " << in.region.floatGlbPhi(calo.hwPhi) << ") qual: " << calo.hwEmID << std::endl;
+                  << "(glb phi: " << in.region.floatGlbPhi(calo.hwPhi) << ") qual: " << std::bitset<4>(calo.hwEmID)
+                  << std::endl;
     }
   }
 
@@ -222,7 +225,15 @@ void PFTkEGAlgoEmulator::eg_algo(const PFRegionEmu &region,
     // check if brem recovery is on
     if (!cfg.doBremRecovery || cfg.writeBeforeBremRecovery) {
       // 1. create EG objects before brem recovery
-      addEgObjsToPF(egstas, egobjs, egeleobjs, emcalo, track, ic, calo.hwEmID, calo.hwPt, itk);
+      unsigned int egQual = calo.hwEmID;
+      // If we write both objects with and without brem-recovery
+      // bit 3 is used for the brem-recovery bit: if set = no recovery
+      // (for consistency with the barrel hwQual where by default the brem recovery is done upstream)
+      if (cfg.writeBeforeBremRecovery && cfg.doBremRecovery) {
+        egQual = calo.hwEmID | 0x8;
+      }
+
+      addEgObjsToPF(egstas, egobjs, egeleobjs, emcalo, track, ic, egQual, calo.hwPt, itk);
     }
 
     if (!cfg.doBremRecovery)
@@ -245,13 +256,13 @@ void PFTkEGAlgoEmulator::eg_algo(const PFRegionEmu &region,
 
     // 2. create EG objects with brem recovery
     // NOTE: duplicating the object is suboptimal but this is done for keeping things as in TDR code...
-    addEgObjsToPF(egstas, egobjs, egeleobjs, emcalo, track, ic, calo.hwEmID + 2, ptBremReco, itk, components);
+    addEgObjsToPF(egstas, egobjs, egeleobjs, emcalo, track, ic, calo.hwEmID, ptBremReco, itk, components);
   }
 }
 
 EGObjEmu &PFTkEGAlgoEmulator::addEGStaToPF(std::vector<EGObjEmu> &egobjs,
                                            const EmCaloObjEmu &calo,
-                                           const int hwQual,
+                                           const unsigned int hwQual,
                                            const pt_t ptCorr,
                                            const std::vector<unsigned int> &components) const {
   EGObjEmu egsta;
@@ -261,40 +272,57 @@ EGObjEmu &PFTkEGAlgoEmulator::addEGStaToPF(std::vector<EGObjEmu> &egobjs,
   egsta.hwPhi = calo.hwPhi;
   egsta.hwQual = hwQual;
   egobjs.push_back(egsta);
+
+  if (debug_ > 2)
+    dbgCout() << "[REF] EGSta pt: " << egsta.hwPt << " eta: " << egsta.hwEta << " phi: " << egsta.hwPhi
+              << " qual: " << std::bitset<4>(egsta.hwQual) << " packed: " << egsta.pack().to_string(16) << std::endl;
+
   return egobjs.back();
 }
 
 EGIsoObjEmu &PFTkEGAlgoEmulator::addEGIsoToPF(std::vector<EGIsoObjEmu> &egobjs,
                                               const EmCaloObjEmu &calo,
-                                              const int hwQual,
+                                              const unsigned int hwQual,
                                               const pt_t ptCorr) const {
   EGIsoObjEmu egiso;
   egiso.clear();
   egiso.hwPt = ptCorr;
   egiso.hwEta = calo.hwEta;
   egiso.hwPhi = calo.hwPhi;
-  egiso.hwQual = hwQual;
+  unsigned int egHwQual = hwQual;
+  if (cfg.doEndcapHwQual) {
+    // 1. zero-suppress the loose EG-ID (bit 1)
+    // 2. for now use the standalone tight definition (bit 0) to set the tight point for photons (bit 2)
+    egHwQual = (hwQual & 0x9) | ((hwQual & 0x1) << 2);
+  }
+  egiso.hwQual = egHwQual;
   egiso.srcCluster = calo.src;
   egobjs.push_back(egiso);
 
   if (debug_ > 2)
     dbgCout() << "[REF] EGIsoObjEmu pt: " << egiso.hwPt << " eta: " << egiso.hwEta << " phi: " << egiso.hwPhi
-              << " qual: " << egiso.hwQual << " packed: " << egiso.pack().to_string(16) << std::endl;
+              << " qual: " << std::bitset<4>(egiso.hwQual) << " packed: " << egiso.pack().to_string(16) << std::endl;
 
   return egobjs.back();
 }
 
 EGIsoEleObjEmu &PFTkEGAlgoEmulator::addEGIsoEleToPF(std::vector<EGIsoEleObjEmu> &egobjs,
                                                     const EmCaloObjEmu &calo,
                                                     const TkObjEmu &track,
-                                                    const int hwQual,
+                                                    const unsigned int hwQual,
                                                     const pt_t ptCorr) const {
   EGIsoEleObjEmu egiso;
   egiso.clear();
   egiso.hwPt = ptCorr;
   egiso.hwEta = calo.hwEta;
   egiso.hwPhi = calo.hwPhi;
-  egiso.hwQual = hwQual;
+  unsigned int egHwQual = hwQual;
+  if (cfg.doEndcapHwQual) {
+    // 1. zero-suppress the loose EG-ID (bit 1)
+    // 2. for now use the standalone tight definition (bit 0) to set the tight point for eles (bit 1)
+    egHwQual = (hwQual & 0x9) | ((hwQual & 0x1) << 1);
+  }
+  egiso.hwQual = egHwQual;
   egiso.hwDEta = track.hwVtxEta() - egiso.hwEta;
   egiso.hwDPhi = abs(track.hwVtxPhi() - egiso.hwPhi);
   egiso.hwZ0 = track.hwZ0;
@@ -305,7 +333,7 @@ EGIsoEleObjEmu &PFTkEGAlgoEmulator::addEGIsoEleToPF(std::vector<EGIsoEleObjEmu>
 
   if (debug_ > 2)
     dbgCout() << "[REF] EGIsoEleObjEmu pt: " << egiso.hwPt << " eta: " << egiso.hwEta << " phi: " << egiso.hwPhi
-              << " qual: " << egiso.hwQual << " packed: " << egiso.pack().to_string(16) << std::endl;
+              << " qual: " << std::bitset<4>(egiso.hwQual) << " packed: " << egiso.pack().to_string(16) << std::endl;
 
   return egobjs.back();
 }
@@ -316,7 +344,7 @@ void PFTkEGAlgoEmulator::addEgObjsToPF(std::vector<EGObjEmu> &egstas,
                                        const std::vector<EmCaloObjEmu> &emcalo,
                                        const std::vector<TkObjEmu> &track,
                                        const int calo_idx,
-                                       const int hwQual,
+                                       const unsigned int hwQual,
                                        const pt_t ptCorr,
                                        const int tk_idx,
                                        const std::vector<unsigned int> &components) const {