HGCal EGamma ID helpers and Phase-II MiniAOD 93x

PhysicsTools/PatAlgos/plugins/PATElectronProducer.cc

@@ -77,6 +77,7 @@ PATElectronProducer::PATElectronProducer(const edm::ParameterSet & iConfig) :
   pfCandidateMultiMapToken_(usePfCandidateMultiMap_ ? consumes<edm::ValueMap<std::vector<reco::PFCandidateRef>>>(iConfig.getParameter<edm::InputTag>( "pfCandidateMultiMap" )) : edm::EDGetTokenT<edm::ValueMap<std::vector<reco::PFCandidateRef>>>()),
   embedPFCandidate_(iConfig.getParameter<bool>( "embedPFCandidate" )),
   // mva input variables
+  addMVAVariables_(iConfig.getParameter<bool>("addMVAVariables")),
   reducedBarrelRecHitCollection_(iConfig.getParameter<edm::InputTag>("reducedBarrelRecHitCollection")),
   reducedBarrelRecHitCollectionToken_(mayConsume<EcalRecHitCollection>(reducedBarrelRecHitCollection_)),
   reducedEndcapRecHitCollection_(iConfig.getParameter<edm::InputTag>("reducedEndcapRecHitCollection")),
@@ -433,9 +434,11 @@ void PATElectronProducer::produce(edm::Event & iEvent, const edm::EventSetup & i
 	    anElectron.setElectronIDs(ids);
 	  }
 
-          // add missing mva variables
-          std::vector<float> vCov = lazyTools.localCovariances(*( itElectron->superCluster()->seed()));
-          anElectron.setMvaVariables(vCov[1], ip3d);
+          if (addMVAVariables_) {
+            // add missing mva variables
+            std::vector<float> vCov = lazyTools.localCovariances(*( itElectron->superCluster()->seed()));
+            anElectron.setMvaVariables(vCov[1], ip3d);
+          }
 	  // PFClusterIso
 	  if (addPFClusterIso_) {
 	    // Get PFCluster Isolation
@@ -667,9 +670,11 @@ void PATElectronProducer::produce(edm::Event & iEvent, const edm::EventSetup & i
 	}
       }
 
-      // add mva variables
-      std::vector<float> vCov = lazyTools.localCovariances(*( itElectron->superCluster()->seed()));
-      anElectron.setMvaVariables(vCov[1], ip3d);
+      if (addMVAVariables_) {
+        // add mva variables
+        std::vector<float> vCov = lazyTools.localCovariances(*( itElectron->superCluster()->seed()));
+        anElectron.setMvaVariables(vCov[1], ip3d);
+      }
       // PFCluster Isolation
       if (addPFClusterIso_) {
 	// Get PFCluster Isolation
@@ -1110,7 +1115,7 @@ void PATElectronProducer::fillDescriptions(edm::ConfigurationDescriptions & desc
 
 
   // electron shapes
-  iDesc.add<bool>("addElectronShapes", true);
+  iDesc.add<bool>("addMVAVariables", true)->setComment("embed extra variables in pat::Electron : sip3d, sigmaIEtaIPhi");
   iDesc.add<edm::InputTag>("reducedBarrelRecHitCollection", edm::InputTag("reducedEcalRecHitsEB"));
   iDesc.add<edm::InputTag>("reducedEndcapRecHitCollection", edm::InputTag("reducedEcalRecHitsEE"));
 

PhysicsTools/PatAlgos/plugins/PATElectronProducer.h

@@ -99,6 +99,7 @@ namespace pat {
       const bool          embedPFCandidate_;
 
       /// mva input variables
+      const bool addMVAVariables_;
       const edm::InputTag reducedBarrelRecHitCollection_;
       const edm::EDGetTokenT<EcalRecHitCollection> reducedBarrelRecHitCollectionToken_;
       const edm::InputTag reducedEndcapRecHitCollection_;

PhysicsTools/PatAlgos/plugins/PATPhotonProducer.cc

@@ -341,14 +341,13 @@ void PATPhotonProducer::produce(edm::Event & iEvent, const edm::EventSetup & iSe
     // set seed energy
     aPhoton.setSeedEnergy( photonRef->superCluster()->seed()->energy() );
 
+    // set input variables for regression energy correction
+    if (saveRegressionData_) {
     EcalRegressionData ecalRegData;
     ecalRegData.fill(*(photonRef->superCluster()),
 		     recHitsEBHandle.product(),recHitsEEHandle.product(),
 		     ecalGeometry_,ecalTopology_,-1);
 
-
-    // set input variables for regression energy correction
-    if (saveRegressionData_) {
     aPhoton.setEMax( ecalRegData.eMax() );
     aPhoton.setE2nd( ecalRegData.e2nd() );
     aPhoton.setE3x3( ecalRegData.e3x3() );

PhysicsTools/PatAlgos/python/producersLayer1/electronProducer_cfi.py

@@ -11,6 +11,7 @@
     usePfCandidateMultiMap = cms.bool( False ),
 
     # collections for mva input variables
+    addMVAVariables = cms.bool( True ),
     reducedBarrelRecHitCollection = cms.InputTag("reducedEcalRecHitsEB"),
     reducedEndcapRecHitCollection = cms.InputTag("reducedEcalRecHitsEE"),
 

PhysicsTools/PatAlgos/python/slimming/MicroEventContent_cff.py

@@ -123,3 +123,8 @@
 cms.untracked.PSet(branch = cms.untracked.string("patJets_slimmedJetsPuppi__*"),splitLevel=cms.untracked.int32(99)),
 cms.untracked.PSet(branch = cms.untracked.string("EcalRecHitsSorted_reducedEgamma_reducedESRecHits_*"),splitLevel=cms.untracked.int32(99)),
 ])
+
+
+_phase2_hgc_extraCommands = ["keep *_slimmedElectronsFromMultiCl_*_*", "keep *_slimmedPhotonsFromMultiCl_*_*"]
+from Configuration.Eras.Modifier_phase2_hgcal_cff import phase2_hgcal
+phase2_hgcal.toModify(MicroEventContentMC, outputCommands = MicroEventContentMC.outputCommands + _phase2_hgc_extraCommands)

PhysicsTools/PatAlgos/python/slimming/miniAOD_tools.py

@@ -333,6 +333,12 @@ def miniAOD_customizeCommon(process):
     # add DetIdAssociatorRecords to EventSetup (for isolatedTracks)
     process.load("TrackingTools.TrackAssociator.DetIdAssociatorESProducer_cff")
 
+    # EGamma objects from HGCal are not yet in GED
+    # so add companion collections for Phase-II MiniAOD production
+    from Configuration.Eras.Modifier_phase2_hgcal_cff import phase2_hgcal
+    process.load("RecoEgamma.EgammaTools.slimmedEgammaFromMultiCl_cff")
+    phase2_hgcal.toModify(task, func=lambda t: t.add(process.slimmedEgammaFromMultiClTask))
+
 
 def miniAOD_customizeMC(process):
     task = getPatAlgosToolsTask(process)

RecoEcal/Configuration/python/RecoEcal_EventContent_cff.py

@@ -84,9 +84,11 @@
         )
 )
 
-_phase2_hgcal_scCommands = ['keep *_particleFlowSuperClusterHGCal_*_*']
+_phase2_hgcal_scCommands = ['keep *_particleFlowSuperClusterHGCal_*_*', 'keep *_particleFlowSuperClusterHGCalFromMultiCl_*_*']
 _phase2_hgcal_scCommandsAOD = ['keep recoSuperClusters_particleFlowSuperClusterHGCal__*',
-                               'keep recoCaloClusters_particleFlowSuperClusterHGCal__*']
+                               'keep recoCaloClusters_particleFlowSuperClusterHGCal__*',
+                               'keep recoSuperClusters_particleFlowSuperClusterHGCalFromMultiCl__*',
+                               'keep recoCaloClusters_particleFlowSuperClusterHGCalFromMultiCl__*']
 _phase2_hgcal_RecoEcalFEVT = RecoEcalFEVT.clone()
 _phase2_hgcal_RecoEcalFEVT.outputCommands += _phase2_hgcal_scCommands
 _phase2_hgcal_RecoEcalRECO = RecoEcalRECO.clone()

RecoEgamma/EgammaTools/BuildFile.xml

@@ -3,6 +3,8 @@
 <use   name="DataFormats/EgammaReco"/>
 <use   name="DataFormats/EgammaCandidates"/>
 <use   name="Geometry/CaloGeometry"/>
+<use   name="FastSimulation/CaloGeometryTools"/>
+<use   name="RecoLocalCalo/HGCalRecAlgos"/>
 <use   name="clhep"/>
 <use   name="root"/>
 <use   name="PhysicsTools/MVAComputer"/>

RecoEgamma/EgammaTools/interface/EgammaPCAHelper.h

@@ -0,0 +1,118 @@
+//--------------------------------------------------------------------------------------------------
+//
+// EGammaPCAHelper
+//
+// Helper Class to compute PCA
+//
+//
+//--------------------------------------------------------------------------------------------------
+#ifndef RecoEgamma_EgammaTools_EGammaPCAHelper_h
+#define RecoEgamma_EgammaTools_EGammaPCAHelper_h
+
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+
+#include "DataFormats/CaloRecHit/interface/CaloCluster.h"
+#include "DataFormats/HGCRecHit/interface/HGCRecHitCollections.h"
+#include "DataFormats/ForwardDetId/interface/HGCalDetId.h"
+#include "DataFormats/Math/interface/Vector3D.h"
+#include "DataFormats/ParticleFlowReco/interface/HGCalMultiCluster.h"
+
+#include "RecoEgamma/EgammaTools/interface/Spot.h"
+#include "RecoEgamma/EgammaTools/interface/LongDeps.h"
+#include "RecoEgamma/EgammaTools/interface/ShowerDepth.h"
+#include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h"
+#include "FastSimulation/CaloGeometryTools/interface/Transform3DPJ.h"
+#include <map>
+
+#include "TPrincipal.h"
+
+class HGCalRecHit;
+
+namespace hgcal {
+
+  class EGammaPCAHelper
+  {
+    public:
+      typedef ROOT::Math::Transform3DPJ Transform3D;
+      typedef ROOT::Math::Transform3DPJ::Point Point;
+
+      EGammaPCAHelper();
+      ~EGammaPCAHelper();
+
+      // for the GsfElectrons
+      void storeRecHits(const reco::CaloCluster & theCluster );
+      void storeRecHits(const reco::HGCalMultiCluster &cluster );
+
+      const TPrincipal & pcaResult();
+      /// to set from outside - once per event
+      void setHitMap( std::map<DetId,const HGCRecHit *> * hitMap) ;
+      /// to compute from inside - once per event
+      void fillHitMap(const HGCRecHitCollection & HGCEERecHits,
+                      const HGCRecHitCollection & HGCFHRecHits,
+                      const HGCRecHitCollection & HGCBHRecHits);
+
+      std::map<DetId,const HGCRecHit *> * getHitMap(){return hitMap_;}
+
+      void setRecHitTools(const hgcal::RecHitTools * recHitTools );
+
+      inline void setdEdXWeights(const std::vector<double> & dEdX){ dEdXWeights_ = dEdX;}
+
+      void pcaInitialComputation() {
+          computePCA(-1.,false);
+      }
+
+      void computePCA(float radius, bool withHalo=true);
+      const math::XYZPoint  & barycenter() const {return barycenter_;}
+      const math::XYZVector & axis() const {return axis_;}
+
+      void computeShowerWidth(float radius, bool withHalo=true);
+
+      inline double sigmaUU() const { return checkIteration()? sigu_ : -1. ;}
+      inline double sigmaVV() const { return checkIteration()? sigv_ : -1. ;}
+      inline double sigmaEE() const { return checkIteration()? sige_ : -1. ;}
+      inline double sigmaPP() const { return checkIteration()? sigp_ : -1. ;}
+
+      inline const TVectorD& eigenValues () const {return *pca_->GetEigenValues();}
+      inline const TVectorD& sigmas() const {return *pca_->GetSigmas();}
+      // contains maxlayer+1 values, first layer is [1]
+      LongDeps  energyPerLayer(float radius, bool withHalo=true) ;
+
+      float clusterDepthCompatibility(const LongDeps &, float & measuredDepth, float& expectedDepth, float& expectedSigma );
+      void printHits( float radius) const;
+      void clear();
+
+    private:
+      bool checkIteration() const ;
+      void storeRecHits(const std::vector<std::pair<DetId, float>> &hf);
+      float findZFirstLayer(const LongDeps&) const;
+
+      bool recHitsStored_;
+      bool debug_;
+
+      //parameters
+      std::vector<double> dEdXWeights_;
+      std::vector<double> invThicknessCorrection_;
+
+      int hitMapOrigin_; // 0 not initialized; 1 set from outside ; 2 set from inside
+      const reco::CaloCluster * theCluster_;
+      std::map<DetId, const HGCRecHit *> * hitMap_;
+      std::vector<Spot> theSpots_;
+      int pcaIteration_;
+
+      // output quantities
+      math::XYZPoint barycenter_;
+      math::XYZVector axis_;
+
+      Transform3D trans_;
+      double sigu_,sigv_,sige_,sigp_;
+
+      // helper
+      std::unique_ptr<TPrincipal> pca_;
+      const hgcal::RecHitTools * recHitTools_;
+      ShowerDepth showerDepth_;
+
+  };
+
+}
+
+#endif

RecoEgamma/EgammaTools/interface/HGCalEgammaIDHelper.h

@@ -0,0 +1,86 @@
+//--------------------------------------------------------------------------------------------------
+//
+// EGammaID Helper
+//
+// Helper Class to compute HGCal Egamma cluster ID observables
+//
+// Authors: F. Beaudette, A. Lobanov, N. Smith
+//--------------------------------------------------------------------------------------------------
+
+#ifndef RecoEgamma_EgammaTools_HGCalEgammaIDHelper_h
+#define RecoEgamma_EgammaTools_HGCalEgammaIDHelper_h
+
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+#include "FWCore/Framework/interface/ESHandle.h"
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Utilities/interface/InputTag.h"
+#include "FWCore/Framework/interface/ConsumesCollector.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+
+
+#include "DataFormats/EgammaCandidates/interface/Photon.h"
+#include "DataFormats/EgammaCandidates/interface/GsfElectronFwd.h"
+#include "DataFormats/CaloRecHit/interface/CaloCluster.h"
+#include "DataFormats/HGCRecHit/interface/HGCRecHitCollections.h"
+
+#include "RecoEgamma/EgammaTools/interface/EgammaPCAHelper.h"
+#include "RecoEgamma/EgammaTools/interface/LongDeps.h"
+#include <vector>
+#include "HGCalIsoCalculator.h"
+
+class HGCalEgammaIDHelper {
+public:
+    HGCalEgammaIDHelper() {}
+    HGCalEgammaIDHelper(const edm::ParameterSet &, edm::ConsumesCollector && iC);
+    ~HGCalEgammaIDHelper() {}
+
+    // Use eventInit once per event
+    void eventInit(const edm::Event& iEvent,const edm::EventSetup &iSetup);
+
+    // Call computeHGCAL before accessing results below
+    void computeHGCAL(const reco::Photon & thePhoton, float radius);
+    void computeHGCAL(const reco::GsfElectron & theElectron, float radius);
+
+    // PCA results
+    double sigmaUU() const {  return pcaHelper_.sigmaUU();}
+    double sigmaVV() const {  return pcaHelper_.sigmaVV();}
+    double sigmaEE() const {  return pcaHelper_.sigmaEE();}
+    double sigmaPP() const {  return pcaHelper_.sigmaPP();}
+    const TVectorD& eigenValues () const {return pcaHelper_.eigenValues();}
+    const TVectorD& sigmas() const {return pcaHelper_.sigmas();}
+    const math::XYZPoint  & barycenter() const {return pcaHelper_.barycenter();}
+    const math::XYZVector & axis() const {return pcaHelper_.axis();}
+
+    // longitudinal energy deposits and energy per subdetector as well as layers crossed
+    hgcal::LongDeps energyPerLayer(float radius, bool withHalo=true) {
+        return pcaHelper_.energyPerLayer(radius,withHalo);
+    }
+
+    // shower depth (distance between start and shower max) from ShowerDepth tool
+    float clusterDepthCompatibility(const hgcal::LongDeps & ld, float & measDepth, float & expDepth, float & expSigma)
+        { return pcaHelper_.clusterDepthCompatibility(ld,measDepth,expDepth,expSigma);}
+
+    // projective calo isolation
+    inline float getIsolationRing(unsigned int ring) const { return isoHelper_.getIso(ring); };
+
+
+    // for debugging purposes
+    void printHits(float radius) const { pcaHelper_.printHits(radius); }
+    const hgcal::EGammaPCAHelper * pcaHelper () const {return &pcaHelper_;}
+
+private:
+    edm::InputTag  eeRecHitInputTag_;
+    edm::InputTag  fhRecHitInputTag_;
+    edm::InputTag  bhRecHitInputTag_;
+
+    std::vector<double> dEdXWeights_;
+    hgcal::EGammaPCAHelper pcaHelper_;
+    HGCalIsoCalculator isoHelper_;
+    edm::EDGetTokenT<HGCRecHitCollection> recHitsEE_;
+    edm::EDGetTokenT<HGCRecHitCollection> recHitsFH_;
+    edm::EDGetTokenT<HGCRecHitCollection> recHitsBH_;
+    hgcal::RecHitTools recHitTools_;
+    bool debug_;
+};
+
+#endif