HEEP V7.0 in 90X

PhysicsTools/PatAlgos/python/slimming/miniAOD_tools.py

@@ -186,6 +186,13 @@ def miniAOD_customizeCommon(process):
         cms.InputTag('reducedEgamma','reducedGedGsfElectrons')
     for idmod in electron_ids:
         setupAllVIDIdsInModule(process,idmod,setupVIDElectronSelection,None,False)
+
+    #heepIDVarValueMaps only exists if HEEP V6.1 or HEEP 7.0 ID has already been loaded
+    if hasattr(process,'heepIDVarValueMaps'):
+        process.heepIDVarValueMaps.elesMiniAOD = cms.InputTag('reducedEgamma','reducedGedGsfElectrons')
+        #force HEEP to use miniAOD (otherwise it'll detect the AOD)
+        process.heepIDVarValueMaps.dataFormat = cms.int32(2)
+
 
     #VID Photon IDs
     photon_ids = ['RecoEgamma.PhotonIdentification.Identification.cutBasedPhotonID_Spring15_25ns_V1_cff',

PhysicsTools/SelectorUtils/python/tools/vid_id_tools.py

@@ -96,6 +96,20 @@ def setupVIDElectronSelection(process,cutflow,patProducer=None,addUserData=True)
         idName = cutflow.idName.value()
         addVIDSelectionToPATProducer(patProducer,'egmGsfElectronIDs',idName,addUserData)
 
+    #we know cutflow has a name otherwise an exception would have been thrown in setupVIDSelection
+    #run this for all heep IDs except V60 standard for which it is not needed
+    #it is needed for V61 and V70
+    if (cutflow.idName.value().find("HEEP")!=-1 and
+        cutflow.idName.value()!="heepElectronID-HEEPV60"):
+
+        #not the ideal way but currently the easiest
+        useMiniAOD = process.egmGsfElectronIDs.physicsObjectSrc == cms.InputTag('slimmedElectrons')
+
+        from RecoEgamma.ElectronIdentification.Identification.heepElectronID_tools import addHEEPProducersToSeq
+        addHEEPProducersToSeq(process=process,seq=process.egmGsfElectronIDSequence,
+                              insertIndex=process.egmGsfElectronIDSequence.index(process.egmGsfElectronIDs),
+                              useMiniAOD=useMiniAOD)
+
 ####
 # Muons
 ####

RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h

@@ -191,6 +191,9 @@ class EcalClusterToolsT {
                 static std::vector<float> roundnessBarrelSuperClusters( const reco::SuperCluster &superCluster ,const EcalRecHitCollection &recHits, int weightedPositionMethod = 0, float energyThreshold = 0.0);
                 static std::vector<float> roundnessBarrelSuperClustersUserExtended( const reco::SuperCluster &superCluster ,const EcalRecHitCollection &recHits, int ieta_delta=0, int iphi_delta=0, float energyRHThresh=0.00000, int weightedPositionMethod=0);
                 static std::vector<float> roundnessSelectedBarrelRecHits(const std::vector<std::pair<const EcalRecHit*,float> >&rhVector, int weightedPositionMethod = 0);
+
+                //works out the number of staturated crystals in 5x5
+                static int nrSaturatedCrysIn5x5(const DetId& id,const EcalRecHitCollection* recHits,const CaloTopology *topology);
         private:
                 struct EcalClusterEnergyDeposition
                 { 
@@ -247,7 +250,8 @@ class EcalClusterToolsT {
                 static std::vector<int> getSeedPosition(const std::vector<std::pair<const EcalRecHit*,float> >&RH_ptrs);
                 static float getSumEnergy(const std::vector<std::pair<const EcalRecHit*,float> >&RH_ptrs_fracs);
                 static float computeWeight(float eRH, float energyTotal, int weightedPositionMethod);
-
+
+
 };
 
 // implementation
@@ -1626,6 +1630,31 @@ std::vector<float> EcalClusterToolsT<noZS>::roundnessSelectedBarrelRecHits( cons
     return shapes;
 
 }
+
+
+template<bool noZS>
+int EcalClusterToolsT<noZS>::nrSaturatedCrysIn5x5(const DetId& id,const EcalRecHitCollection* recHits,const CaloTopology *topology)
+{
+  int nrSat=0;
+  CaloNavigator<DetId> cursor = CaloNavigator<DetId>( id, topology->getSubdetectorTopology( id ) );
+
+  for ( int eastNr = -2; eastNr <= 2; ++eastNr ) { //east is eta in barrel
+    for ( int northNr = -2; northNr <= 2; ++northNr ) { //north is phi in barrel
+      cursor.home();
+      cursor.offsetBy( eastNr, northNr);
+      DetId id = *cursor;
+      auto recHitIt = recHits->find(id);
+      if(recHitIt!=recHits->end() && 
+	 recHitIt->checkFlag(EcalRecHit::kSaturated)){
+	nrSat++;
+      }
+
+    }
+  }
+  return nrSat;
+}
+
+
 //private functions useful for roundnessBarrelSuperClusters etc.
 //compute delta iphi between a seed and a particular recHit
 //iphi [1,360]
@@ -1688,6 +1717,7 @@ float EcalClusterToolsT<noZS>::getSumEnergy(const std::vector<std::pair<const Ec
     return sumE;
 }
 
+
 typedef EcalClusterToolsT<false> EcalClusterTools;
 
 namespace noZS {

RecoEgamma/EgammaIsolationAlgos/BuildFile.xml

@@ -8,6 +8,7 @@
 <use   name="DataFormats/RecoCandidate"/>
 <use   name="DataFormats/ParticleFlowReco"/>
 <use   name="DataFormats/ParticleFlowCandidate"/>
+<use   name="DataFormats/PatCandidates"/>
 <use   name="CondFormats/EcalObjects"/>
 <use   name="CondFormats/DataRecord"/>
 <use   name="RecoLocalCalo/EcalRecAlgos"/>

RecoEgamma/EgammaIsolationAlgos/interface/EleTkIsolFromCands.h

@@ -0,0 +1,77 @@
+#ifndef RECOEGAMMA_EGAMMAISOLATIONALGOS_ELETKISOLFROMCANDS_H
+#define RECOEGAMMA_EGAMMAISOLATIONALGOS_ELETKISOLFROMCANDS_H
+
+#include "DataFormats/TrackReco/interface/TrackBase.h"
+#include "DataFormats/PatCandidates/interface/PackedCandidate.h"
+#include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
+#include "DataFormats/Common/interface/View.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
+
+class EleTkIsolFromCands {
+
+private:
+  struct TrkCuts {
+    float minPt;
+    float minDR2;
+    float maxDR2;
+    float minDEta;
+    float maxDZ;
+    float minHits;
+    float minPixelHits;
+    float maxDPtPt;
+    std::vector<reco::TrackBase::TrackQuality> allowedQualities;
+    std::vector<reco::TrackBase::TrackAlgorithm> algosToReject;
+    explicit TrkCuts(const edm::ParameterSet& para);
+    static edm::ParameterSetDescription pSetDescript();
+  };
+
+  TrkCuts barrelCuts_,endcapCuts_;
+
+public:
+  explicit EleTkIsolFromCands(const edm::ParameterSet& para);
+  EleTkIsolFromCands(const EleTkIsolFromCands&)=default;
+  ~EleTkIsolFromCands()=default;
+  EleTkIsolFromCands& operator=(const EleTkIsolFromCands&)=default;
+
+  static edm::ParameterSetDescription pSetDescript();
+
+  std::pair<int,double> calIsol(const reco::TrackBase& trk,const pat::PackedCandidateCollection& cands,const edm::View<reco::GsfElectron>& eles);
+
+  std::pair<int,double> calIsol(const double eleEta,const double elePhi,const double eleVZ,
+				const pat::PackedCandidateCollection& cands,
+				const edm::View<reco::GsfElectron>& eles);
+
+  double calIsolPt(const reco::TrackBase& trk,const pat::PackedCandidateCollection& cands,
+		   const edm::View<reco::GsfElectron>& eles){
+    return calIsol(trk,cands,eles).second;
+  }
+
+  double calIsolPt(const double eleEta,const double elePhi,const double eleVZ,
+		   const pat::PackedCandidateCollection& cands,
+		   const edm::View<reco::GsfElectron>& eles){
+    return calIsol(eleEta,elePhi,eleVZ,cands,eles).second;
+  }
+
+  static bool passTrkSel(const reco::Track& trk,
+			 const double trkPt,
+			 const TrkCuts& cuts,
+			 const double eleEta,const double elePhi,
+			 const double eleVZ);
+
+private:
+  //no qualities specified, accept all, ORed
+  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);
+  //for PF electron candidates the "trk pt" is not the track pt
+  //its the trk-calo gsfele combination energy * trk sin(theta)
+  //so the idea is to match with the gsf electron and get the orginal
+  //gsftrack's pt
+  double getTrkPt(const reco::TrackBase& trk,
+		  const edm::View<reco::GsfElectron>& eles);
+};
+
+
+#endif

RecoEgamma/EgammaIsolationAlgos/src/EleTkIsolFromCands.cc

@@ -0,0 +1,161 @@
+#include "RecoEgamma/EgammaIsolationAlgos/interface/EleTkIsolFromCands.h"
+
+#include "DataFormats/Math/interface/deltaR.h"
+#include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
+
+EleTkIsolFromCands::TrkCuts::TrkCuts(const edm::ParameterSet& para)
+{
+  minPt = para.getParameter<double>("minPt");
+  auto sq = [](double val){return val*val;};
+  minDR2 = sq(para.getParameter<double>("minDR"));
+  maxDR2 = sq(para.getParameter<double>("maxDR"));
+  minDEta = para.getParameter<double>("minDEta");
+  maxDZ = para.getParameter<double>("maxDZ");
+  minHits = para.getParameter<int>("minHits");
+  minPixelHits = para.getParameter<int>("minPixelHits");
+  maxDPtPt = para.getParameter<double>("maxDPtPt");
+
+  auto qualNames = para.getParameter<std::vector<std::string> >("allowedQualities");
+  auto algoNames = para.getParameter<std::vector<std::string> >("algosToReject");
+
+  for(auto& qualName : qualNames){
+    allowedQualities.push_back(reco::TrackBase::qualityByName(qualName));
+  }
+  for(auto& algoName : algoNames){
+    algosToReject.push_back(reco::TrackBase::algoByName(algoName));
+  }
+  std::sort(algosToReject.begin(),algosToReject.end());
+
+}
+
+edm::ParameterSetDescription EleTkIsolFromCands::TrkCuts::pSetDescript()
+{
+  edm::ParameterSetDescription desc;
+  desc.add<double>("minPt",1.0);
+  desc.add<double>("maxDR",0.3);
+  desc.add<double>("minDR",0.000);
+  desc.add<double>("minDEta",0.005);
+  desc.add<double>("maxDZ",0.1);
+  desc.add<double>("maxDPtPt",-1);
+  desc.add<int>("minHits",8);
+  desc.add<int>("minPixelHits",1);
+  desc.add<std::vector<std::string> >("allowedQualities");
+  desc.add<std::vector<std::string> >("algosToReject");
+  return desc;
+}
+
+EleTkIsolFromCands::EleTkIsolFromCands(const edm::ParameterSet& para):
+  barrelCuts_(para.getParameter<edm::ParameterSet>("barrelCuts")),
+  endcapCuts_(para.getParameter<edm::ParameterSet>("endcapCuts"))          
+{
+
+
+}
+
+edm::ParameterSetDescription EleTkIsolFromCands::pSetDescript()
+{
+  edm::ParameterSetDescription desc;
+  desc.add("barrelCuts",TrkCuts::pSetDescript());
+  desc.add("endcapCuts",TrkCuts::pSetDescript());
+  return desc;
+}
+
+std::pair<int,double> 
+EleTkIsolFromCands::calIsol(const reco::TrackBase& eleTrk,
+			    const pat::PackedCandidateCollection& cands,
+			    const edm::View<reco::GsfElectron>& eles)
+{
+  return calIsol(eleTrk.eta(),eleTrk.phi(),eleTrk.vz(),cands,eles);
+}
+
+std::pair<int,double> 
+EleTkIsolFromCands::calIsol(const double eleEta,const double elePhi,
+			    const double eleVZ,
+			    const pat::PackedCandidateCollection& cands,
+			    const edm::View<reco::GsfElectron>& eles)
+{
+
+  double ptSum=0.;
+  int nrTrks=0;
+
+  const TrkCuts& cuts = std::abs(eleEta)<1.5 ? barrelCuts_ : endcapCuts_;
+
+  for(auto& cand  : cands){
+    if(cand.charge()!=0){
+      const reco::Track& trk = cand.pseudoTrack(); 
+      double trkPt = std::abs(cand.pdgId())!=11 ? trk.pt() : getTrkPt(trk,eles);
+      if(passTrkSel(trk,trkPt,cuts,eleEta,elePhi,eleVZ)){	
+	ptSum+=trkPt;
+	nrTrks++;
+      }
+    }
+  }
+  return {nrTrks,ptSum};	
+}
+
+
+bool EleTkIsolFromCands::passTrkSel(const reco::Track& trk,
+				    const double trkPt,const TrkCuts& cuts,
+				    const double eleEta,const double elePhi,
+				    const double eleVZ)
+{
+  const float dR2 = reco::deltaR2(eleEta,elePhi,trk.eta(),trk.phi());
+  const float dEta = trk.eta()-eleEta;
+  const float dZ = eleVZ - trk.vz();
+
+  return dR2>=cuts.minDR2 && dR2<=cuts.maxDR2 && 
+    std::abs(dEta)>=cuts.minDEta && 
+    std::abs(dZ)<cuts.maxDZ &&
+    trk.hitPattern().numberOfValidHits() >= cuts.minHits &&
+    trk.hitPattern().numberOfValidPixelHits() >=cuts.minPixelHits &&
+    (trk.ptError()/trkPt < cuts.maxDPtPt || cuts.maxDPtPt<0) && 
+    passQual(trk,cuts.allowedQualities) &&
+    passAlgo(trk,cuts.algosToReject) &&
+    trkPt > cuts.minPt;
+}
+
+
+
+bool EleTkIsolFromCands::
+passQual(const reco::TrackBase& trk,
+	 const std::vector<reco::TrackBase::TrackQuality>& quals)
+{
+  if(quals.empty()) return true;
+
+  for(auto qual : quals) {
+    if(trk.quality(qual)) return true;
+  }
+
+  return false;  
+}
+
+bool EleTkIsolFromCands::
+passAlgo(const reco::TrackBase& trk,
+	 const std::vector<reco::TrackBase::TrackAlgorithm>& algosToRej)
+{
+  return algosToRej.empty() || !std::binary_search(algosToRej.begin(),algosToRej.end(),trk.algo());
+}
+
+//so the working theory here is that the track we have is the electrons gsf track
+//if so, lets get the pt of the gsf track before E/p combinations
+//if no match found to a gsf ele with a gsftrack, return the pt of the input track
+double EleTkIsolFromCands::
+getTrkPt(const reco::TrackBase& trk,
+	 const edm::View<reco::GsfElectron>& eles)
+{
+  //note, the trk.eta(),trk.phi() should be identical to the gsf track eta,phi
+  //although this may not be the case due to roundings after packing
+  auto match=[](const reco::TrackBase& trk,const reco::GsfElectron& ele){
+    return std::abs(trk.eta()-ele.gsfTrack()->eta())<0.001 &&
+    reco::deltaPhi(trk.phi(),ele.gsfTrack()->phi())<0.001;// && 
+  };
+  for(auto& ele : eles){
+    if(ele.gsfTrack().isNonnull()){
+      if(match(trk,ele)){
+	return ele.gsfTrack()->pt();
+      }
+    }
+  }
+  return trk.pt();
+
+}

RecoEgamma/ElectronIdentification/BuildFile.xml

@@ -7,6 +7,7 @@
 <use   name="MagneticField/Records"/>
 <use   name="MagneticField/Engine"/>
 <use   name="RecoEgamma/EgammaTools"/>
+<use   name="RecoEgamma/EgammaIsolationAlgos"/>
 <use   name="PhysicsTools/SelectorUtils"/>
 
 <export>

RecoEgamma/ElectronIdentification/interface/EBEECutValues.h

@@ -9,21 +9,27 @@ namespace reco {
   typedef edm::Ptr<reco::GsfElectron> GsfElectronPtr;
 }
 
-class EBEECutValues {
+template<typename T>
+class EBEECutValuesT {
 private:
-  double barrel_;
-  double endcap_;
+  T barrel_;
+  T endcap_;
   const double barrelCutOff_=1.479; //this is currrently used to identify if object is barrel or endcap but may change
 
 public:
-  EBEECutValues(const edm::ParameterSet& params,const std::string& name):
-    barrel_(params.getParameter<double>(name+"EB")),
-    endcap_(params.getParameter<double>(name+"EE")){}
-  double operator()(const reco::GsfElectronPtr& cand)const{return isBarrel(cand) ? barrel_ : endcap_;}
+  EBEECutValuesT(const edm::ParameterSet& params,const std::string& name):
+    barrel_(params.getParameter<T>(name+"EB")),
+    endcap_(params.getParameter<T>(name+"EE")){}
+  T operator()(const reco::GsfElectronPtr& cand)const{return isBarrel(cand) ? barrel_ : endcap_;}
 
 private:
   const bool isBarrel(const reco::GsfElectronPtr& cand)const{return std::abs(cand->superCluster()->position().eta())<barrelCutOff_;}
 
 };
 
+typedef EBEECutValuesT<double> EBEECutValues;
+typedef EBEECutValuesT<int> EBEECutValuesInt;
+
+
+
 #endif