GsfElectronAlgo.cc

Go to the documentation of this file.

#include "DataFormats/EcalDetId/interface/EBDetId.h"
#include "DataFormats/EcalDetId/interface/EEDetId.h"
#include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
#include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
#include "DataFormats/EgammaCandidates/interface/GsfElectronFwd.h"
#include "DataFormats/EgammaReco/interface/BasicCluster.h"
#include "DataFormats/EgammaReco/interface/ElectronSeed.h"
#include "DataFormats/ForwardDetId/interface/ForwardSubdetector.h"
#include "DataFormats/ForwardDetId/interface/HGCalDetId.h"
#include "DataFormats/GeometryVector/interface/GlobalPoint.h"
#include "DataFormats/GeometryVector/interface/GlobalVector.h"
#include "DataFormats/Math/interface/LorentzVector.h"
#include "DataFormats/ParticleFlowReco/interface/GsfPFRecTrack.h"
#include "DataFormats/TrackingRecHit/interface/TrackingRecHitFwd.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/Utilities/interface/isFinite.h"
#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
#include "RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h"
#include "RecoEcal/EgammaCoreTools/interface/EcalTools.h"
#include "RecoEgamma/EgammaElectronAlgos/interface/EgAmbiguityTools.h"
#include "RecoEgamma/EgammaElectronAlgos/interface/ElectronClassification.h"
#include "RecoEgamma/EgammaElectronAlgos/interface/ElectronEnergyCorrector.h"
#include "RecoEgamma/EgammaElectronAlgos/interface/ElectronMomentumCorrector.h"
#include "RecoEgamma/EgammaElectronAlgos/interface/ElectronUtilities.h"
#include "RecoEgamma/EgammaElectronAlgos/interface/GsfElectronAlgo.h"
#include "RecoEgamma/EgammaElectronAlgos/interface/ecalClusterEnergyUncertaintyElectronSpecific.h"
#include "CommonTools/Egamma/interface/ConversionTools.h"
#include "RecoEcal/EgammaCoreTools/interface/EgammaLocalCovParamDefaults.h"
#include "CondFormats/DataRecord/interface/HcalPFCutsRcd.h"
#include "CondTools/Hcal/interface/HcalPFCutsHandler.h"
#include "Geometry/CaloTopology/interface/HcalTopology.h"

#include <Math/Point3D.h>
#include <memory>

#include <algorithm>
#include <sstream>

using namespace edm;
using namespace reco;

GsfElectronAlgo::HeavyObjectCache::HeavyObjectCache(const edm::ParameterSet& conf) {
  // soft electron MVA
  SoftElectronMVAEstimator::Configuration sconfig;
  sconfig.vweightsfiles = conf.getParameter<std::vector<std::string>>("SoftElecMVAFilesString");
  sElectronMVAEstimator = std::make_unique<SoftElectronMVAEstimator>(sconfig);
  // isolated electron MVA
  ElectronMVAEstimator::Configuration iconfig;
  iconfig.vweightsfiles = conf.getParameter<std::vector<std::string>>("ElecMVAFilesString");
  iElectronMVAEstimator = std::make_unique<ElectronMVAEstimator>(iconfig);

  // Here we will have to load the DNN PFID if present in the config
  egammaTools::DNNConfiguration dconfig;
  const auto& pset_dnn = conf.getParameter<edm::ParameterSet>("EleDNNPFid");
  const bool dnnEnabled = pset_dnn.getParameter<bool>("enabled");
  if (dnnEnabled) {
    dconfig.inputTensorName = pset_dnn.getParameter<std::string>("inputTensorName");
    dconfig.outputTensorName = pset_dnn.getParameter<std::string>("outputTensorName");
    dconfig.modelsFiles = pset_dnn.getParameter<std::vector<std::string>>("modelsFiles");
    dconfig.scalersFiles = pset_dnn.getParameter<std::vector<std::string>>("scalersFiles");
    dconfig.outputDim = pset_dnn.getParameter<std::vector<unsigned int>>("outputDim");
    const auto useEBModelInGap = pset_dnn.getParameter<bool>("useEBModelInGap");
    iElectronDNNEstimator = std::make_unique<ElectronDNNEstimator>(dconfig, useEBModelInGap);
  }
}

//===================================================================
// GsfElectronAlgo::EventData
//===================================================================

struct GsfElectronAlgo::EventData {
  // utilities
  void retreiveOriginalTrackCollections(const reco::TrackRef&, const reco::GsfTrackRef&);

  // general
  edm::Event const* event;
  const reco::BeamSpot* beamspot;

  // input collections
  edm::Handle<reco::GsfElectronCoreCollection> coreElectrons;
  edm::Handle<EcalRecHitCollection> barrelRecHits;
  edm::Handle<EcalRecHitCollection> endcapRecHits;
  edm::Handle<reco::TrackCollection> currentCtfTracks;
  edm::Handle<reco::ElectronSeedCollection> seeds;
  edm::Handle<reco::VertexCollection> vertices;
  edm::Handle<reco::ConversionCollection> conversions;

  // isolation helpers
  EgammaHcalIsolation hadIsolation03, hadIsolation04;
  EgammaHcalIsolation hadIsolation03Bc, hadIsolation04Bc;
  EgammaRecHitIsolation ecalBarrelIsol03, ecalBarrelIsol04;
  EgammaRecHitIsolation ecalEndcapIsol03, ecalEndcapIsol04;

  EleTkIsolFromCands tkIsol03Calc;
  EleTkIsolFromCands tkIsol04Calc;
  EleTkIsolFromCands tkIsolHEEP03Calc;
  EleTkIsolFromCands tkIsolHEEP04Calc;

  edm::Handle<reco::TrackCollection> originalCtfTracks;
  edm::Handle<reco::GsfTrackCollection> originalGsfTracks;

  bool originalCtfTrackCollectionRetreived = false;
  bool originalGsfTrackCollectionRetreived = false;

  //PFCluster Isolation handles
  edm::Handle<reco::PFClusterCollection> ecalClustersHandle;
  std::vector<edm::Handle<reco::PFClusterCollection>> hcalClustersHandle;
};

//===================================================================
// GsfElectronAlgo::ElectronData
//===================================================================

struct GsfElectronAlgo::ElectronData {
  // Refs to subproducts
  const reco::GsfElectronCoreRef coreRef;
  const reco::GsfTrackRef gsfTrackRef;
  const reco::SuperClusterRef superClusterRef;
  reco::TrackRef ctfTrackRef;
  float shFracInnerHits;
  const reco::BeamSpot beamSpot;

  // constructors
  ElectronData(const reco::GsfElectronCoreRef& core, const reco::BeamSpot& bs);

  // utilities
  void computeCharge(int& charge, reco::GsfElectron::ChargeInfo& info);
  reco::CaloClusterPtr getEleBasicCluster(MultiTrajectoryStateTransform const&);
  bool calculateTSOS(MultiTrajectoryStateTransform const&, GsfConstraintAtVertex const&);
  void calculateMode();
  reco::Candidate::LorentzVector calculateMomentum();

  // TSOS
  TrajectoryStateOnSurface innTSOS;
  TrajectoryStateOnSurface outTSOS;
  TrajectoryStateOnSurface vtxTSOS;
  TrajectoryStateOnSurface sclTSOS;
  TrajectoryStateOnSurface seedTSOS;
  TrajectoryStateOnSurface eleTSOS;
  TrajectoryStateOnSurface constrainedVtxTSOS;

  // mode
  GlobalVector innMom, seedMom, eleMom, sclMom, vtxMom, outMom;
  GlobalPoint innPos, seedPos, elePos, sclPos, vtxPos, outPos;
  GlobalVector vtxMomWithConstraint;
};

void GsfElectronAlgo::EventData::retreiveOriginalTrackCollections(const reco::TrackRef& ctfTrack,
                                                                  const reco::GsfTrackRef& gsfTrack) {
  if ((!originalCtfTrackCollectionRetreived) && (ctfTrack.isNonnull())) {
    event->get(ctfTrack.id(), originalCtfTracks);
    originalCtfTrackCollectionRetreived = true;
  }
  if ((!originalGsfTrackCollectionRetreived) && (gsfTrack.isNonnull())) {
    event->get(gsfTrack.id(), originalGsfTracks);
    originalGsfTrackCollectionRetreived = true;
  }
}

GsfElectronAlgo::ElectronData::ElectronData(const reco::GsfElectronCoreRef& core, const reco::BeamSpot& bs)
    : coreRef(core),
      gsfTrackRef(coreRef->gsfTrack()),
      superClusterRef(coreRef->superCluster()),
      ctfTrackRef(coreRef->ctfTrack()),
      shFracInnerHits(coreRef->ctfGsfOverlap()),
      beamSpot(bs) {}

void GsfElectronAlgo::ElectronData::computeCharge(int& charge, GsfElectron::ChargeInfo& info) {
  // determine charge from SC
  GlobalPoint orig, scpos;
  ele_convert(beamSpot.position(), orig);
  ele_convert(superClusterRef->position(), scpos);
  GlobalVector scvect(scpos - orig);
  GlobalPoint inntkpos = innTSOS.globalPosition();
  GlobalVector inntkvect = GlobalVector(inntkpos - orig);
  float dPhiInnEle = normalizedPhi(scvect.barePhi() - inntkvect.barePhi());
  if (dPhiInnEle > 0)
    info.scPixCharge = -1;
  else
    info.scPixCharge = 1;

  // flags
  int chargeGsf = gsfTrackRef->charge();
  info.isGsfScPixConsistent = ((chargeGsf * info.scPixCharge) > 0);
  info.isGsfCtfConsistent = (ctfTrackRef.isNonnull() && ((chargeGsf * ctfTrackRef->charge()) > 0));
  info.isGsfCtfScPixConsistent = (info.isGsfScPixConsistent && info.isGsfCtfConsistent);

  // default charge
  if (info.isGsfScPixConsistent || ctfTrackRef.isNull()) {
    charge = info.scPixCharge;
  } else {
    charge = ctfTrackRef->charge();
  }
}

CaloClusterPtr GsfElectronAlgo::ElectronData::getEleBasicCluster(MultiTrajectoryStateTransform const& mtsTransform) {
  CaloClusterPtr eleRef;
  TrajectoryStateOnSurface tempTSOS;
  TrajectoryStateOnSurface outTSOS = mtsTransform.outerStateOnSurface(*gsfTrackRef);
  float dphimin = 1.e30;
  for (auto const& bc : superClusterRef->clusters()) {
    GlobalPoint posclu(bc->position().x(), bc->position().y(), bc->position().z());
    tempTSOS = mtsTransform.extrapolatedState(outTSOS, posclu);
    if (!tempTSOS.isValid())
      tempTSOS = outTSOS;
    GlobalPoint extrap = tempTSOS.globalPosition();
    float dphi = EleRelPointPair(posclu, extrap, beamSpot.position()).dPhi();
    if (std::abs(dphi) < dphimin) {
      dphimin = std::abs(dphi);
      eleRef = bc;
      eleTSOS = tempTSOS;
    }
  }
  return eleRef;
}

bool GsfElectronAlgo::ElectronData::calculateTSOS(MultiTrajectoryStateTransform const& mtsTransform,
                                                  GsfConstraintAtVertex const& constraintAtVtx) {
  //at innermost point
  innTSOS = mtsTransform.innerStateOnSurface(*gsfTrackRef);
  if (!innTSOS.isValid())
    return false;

  //at vertex
  // innermost state propagation to the beam spot position
  GlobalPoint bsPos;
  ele_convert(beamSpot.position(), bsPos);
  vtxTSOS = mtsTransform.extrapolatedState(innTSOS, bsPos);
  if (!vtxTSOS.isValid())
    vtxTSOS = innTSOS;

  //at seed
  outTSOS = mtsTransform.outerStateOnSurface(*gsfTrackRef);
  if (!outTSOS.isValid())
    return false;

  //    TrajectoryStateOnSurface seedTSOS
  seedTSOS = mtsTransform.extrapolatedState(outTSOS,
                                            GlobalPoint(superClusterRef->seed()->position().x(),
                                                        superClusterRef->seed()->position().y(),
                                                        superClusterRef->seed()->position().z()));
  if (!seedTSOS.isValid())
    seedTSOS = outTSOS;

  // at scl
  sclTSOS = mtsTransform.extrapolatedState(
      innTSOS, GlobalPoint(superClusterRef->x(), superClusterRef->y(), superClusterRef->z()));
  if (!sclTSOS.isValid())
    sclTSOS = outTSOS;

  // constrained momentum
  constrainedVtxTSOS = constraintAtVtx.constrainAtBeamSpot(*gsfTrackRef, beamSpot);

  return true;
}

void GsfElectronAlgo::ElectronData::calculateMode() {
  multiTrajectoryStateMode::momentumFromModeCartesian(innTSOS, innMom);
  multiTrajectoryStateMode::positionFromModeCartesian(innTSOS, innPos);
  multiTrajectoryStateMode::momentumFromModeCartesian(seedTSOS, seedMom);
  multiTrajectoryStateMode::positionFromModeCartesian(seedTSOS, seedPos);
  multiTrajectoryStateMode::momentumFromModeCartesian(eleTSOS, eleMom);
  multiTrajectoryStateMode::positionFromModeCartesian(eleTSOS, elePos);
  multiTrajectoryStateMode::momentumFromModeCartesian(sclTSOS, sclMom);
  multiTrajectoryStateMode::positionFromModeCartesian(sclTSOS, sclPos);
  multiTrajectoryStateMode::momentumFromModeCartesian(vtxTSOS, vtxMom);
  multiTrajectoryStateMode::positionFromModeCartesian(vtxTSOS, vtxPos);
  multiTrajectoryStateMode::momentumFromModeCartesian(outTSOS, outMom);
  multiTrajectoryStateMode::positionFromModeCartesian(outTSOS, outPos);
  multiTrajectoryStateMode::momentumFromModeCartesian(constrainedVtxTSOS, vtxMomWithConstraint);
}

Candidate::LorentzVector GsfElectronAlgo::ElectronData::calculateMomentum() {
  double scale = superClusterRef->energy() / vtxMom.mag();
  return Candidate::LorentzVector(
      vtxMom.x() * scale, vtxMom.y() * scale, vtxMom.z() * scale, superClusterRef->energy());
}

reco::GsfElectron::SaturationInfo GsfElectronAlgo::calculateSaturationInfo(const reco::SuperClusterRef& theClus,
                                                                           EventData const& eventData) const {
  reco::GsfElectron::SaturationInfo si;

  const reco::CaloCluster& seedCluster = *(theClus->seed());
  DetId seedXtalId = seedCluster.seed();
  int detector = seedXtalId.subdetId();
  const EcalRecHitCollection* ecalRecHits = nullptr;
  if (detector == EcalBarrel)
    ecalRecHits = eventData.barrelRecHits.product();
  else
    ecalRecHits = eventData.endcapRecHits.product();

  int nSaturatedXtals = 0;
  bool isSeedSaturated = false;
  for (auto&& hitFractionPair : theClus->hitsAndFractions()) {
    auto&& ecalRecHit = ecalRecHits->find(hitFractionPair.first);
    if (ecalRecHit == ecalRecHits->end())
      continue;
    if (ecalRecHit->checkFlag(EcalRecHit::Flags::kSaturated)) {
      nSaturatedXtals++;
      if (seedXtalId == ecalRecHit->detid())
        isSeedSaturated = true;
    }
  }
  si.nSaturatedXtals = nSaturatedXtals;
  si.isSeedSaturated = isSeedSaturated;

  return si;
}

template <bool full5x5>
reco::GsfElectron::ShowerShape GsfElectronAlgo::calculateShowerShape(const reco::SuperClusterRef& theClus,
                                                                     ElectronHcalHelper const& hcalHelperCone,
                                                                     ElectronHcalHelper const& hcalHelperBc,
                                                                     EventData const& eventData,
                                                                     CaloTopology const& topology,
                                                                     CaloGeometry const& geometry,
                                                                     EcalPFRecHitThresholds const& thresholds,
                                                                     const HcalPFCuts* hcalCuts) const {
  using ClusterTools = EcalClusterToolsT<full5x5>;
  reco::GsfElectron::ShowerShape showerShape;

  const reco::CaloCluster& seedCluster = *(theClus->seed());
  // temporary, till CaloCluster->seed() is made available
  DetId seedXtalId = seedCluster.hitsAndFractions()[0].first;
  int detector = seedXtalId.subdetId();

  const EcalRecHitCollection* recHits = nullptr;
  std::vector<int> recHitFlagsToBeExcluded;
  std::vector<int> recHitSeverityToBeExcluded;
  if (detector == EcalBarrel) {
    recHits = eventData.barrelRecHits.product();
    recHitFlagsToBeExcluded = cfg_.recHits.recHitFlagsToBeExcludedBarrel;
    recHitSeverityToBeExcluded = cfg_.recHits.recHitSeverityToBeExcludedBarrel;
  } else {
    recHits = eventData.endcapRecHits.product();
    recHitFlagsToBeExcluded = cfg_.recHits.recHitFlagsToBeExcludedEndcaps;
    recHitSeverityToBeExcluded = cfg_.recHits.recHitSeverityToBeExcludedEndcaps;
  }

  const auto& covariances = ClusterTools::covariances(seedCluster, recHits, &topology, &geometry);

  // do noise-cleaning for full5x5, by passing per crystal PF recHit thresholds and mult values
  // mult values for EB and EE were obtained by dedicated studies
  const auto& localCovariances = full5x5 ? ClusterTools::localCovariances(seedCluster,
                                                                          recHits,
                                                                          &topology,
                                                                          EgammaLocalCovParamDefaults::kRelEnCut,
                                                                          &thresholds,
                                                                          cfg_.cuts.multThresEB,
                                                                          cfg_.cuts.multThresEE)
                                         : ClusterTools::localCovariances(seedCluster, recHits, &topology);

  showerShape.sigmaEtaEta = sqrt(covariances[0]);
  showerShape.sigmaIetaIeta = sqrt(localCovariances[0]);
  if (!edm::isNotFinite(localCovariances[2]))
    showerShape.sigmaIphiIphi = sqrt(localCovariances[2]);
  showerShape.e1x5 = ClusterTools::e1x5(seedCluster, recHits, &topology);
  showerShape.e2x5Max = ClusterTools::e2x5Max(seedCluster, recHits, &topology);
  showerShape.e5x5 = ClusterTools::e5x5(seedCluster, recHits, &topology);
  showerShape.r9 = ClusterTools::e3x3(seedCluster, recHits, &topology) / theClus->rawEnergy();

  const float scale = full5x5 ? showerShape.e5x5 : theClus->energy();

  for (uint id = 0; id < showerShape.hcalOverEcal.size(); ++id) {
    showerShape.hcalOverEcal[id] = hcalHelperCone.hcalESum(*theClus, id + 1, hcalCuts) / scale;
    showerShape.hcalOverEcalBc[id] = hcalHelperBc.hcalESum(*theClus, id + 1, hcalCuts) / scale;
  }
  showerShape.invalidHcal = !hcalHelperBc.hasActiveHcal(*theClus);
  showerShape.hcalTowersBehindClusters = hcalHelperBc.hcalTowersBehindClusters(*theClus);
  showerShape.pre7DepthHcal = false;

  // extra shower shapes
  const float see_by_spp = showerShape.sigmaIetaIeta * showerShape.sigmaIphiIphi;
  if (see_by_spp > 0) {
    showerShape.sigmaIetaIphi = localCovariances[1] / see_by_spp;
  } else if (localCovariances[1] > 0) {
    showerShape.sigmaIetaIphi = 1.f;
  } else {
    showerShape.sigmaIetaIphi = -1.f;
  }
  showerShape.eMax = ClusterTools::eMax(seedCluster, recHits);
  showerShape.e2nd = ClusterTools::e2nd(seedCluster, recHits);
  showerShape.eTop = ClusterTools::eTop(seedCluster, recHits, &topology);
  showerShape.eLeft = ClusterTools::eLeft(seedCluster, recHits, &topology);
  showerShape.eRight = ClusterTools::eRight(seedCluster, recHits, &topology);
  showerShape.eBottom = ClusterTools::eBottom(seedCluster, recHits, &topology);

  showerShape.e2x5Left = ClusterTools::e2x5Left(seedCluster, recHits, &topology);
  showerShape.e2x5Right = ClusterTools::e2x5Right(seedCluster, recHits, &topology);
  showerShape.e2x5Top = ClusterTools::e2x5Top(seedCluster, recHits, &topology);
  showerShape.e2x5Bottom = ClusterTools::e2x5Bottom(seedCluster, recHits, &topology);

  return showerShape;
}

//===================================================================
// GsfElectronAlgo
//===================================================================

GsfElectronAlgo::GsfElectronAlgo(const Tokens& input,
                                 const StrategyConfiguration& strategy,
                                 const CutsConfiguration& cuts,
                                 const ElectronHcalHelper::Configuration& hcalCone,
                                 const ElectronHcalHelper::Configuration& hcalBc,
                                 const IsolationConfiguration& iso,
                                 const PFClusterIsolationConfiguration& pfiso,
                                 const EcalRecHitsConfiguration& recHits,
                                 std::unique_ptr<EcalClusterFunctionBaseClass>&& crackCorrectionFunction,
                                 const RegressionHelper::Configuration& reg,
                                 const edm::ParameterSet& tkIsol03,
                                 const edm::ParameterSet& tkIsol04,
                                 const edm::ParameterSet& tkIsolHEEP03,
                                 const edm::ParameterSet& tkIsolHEEP04,
                                 edm::ConsumesCollector&& cc)
    : cfg_{input, strategy, cuts, iso, pfiso, recHits},
      tkIsol03CalcCfg_(tkIsol03),
      tkIsol04CalcCfg_(tkIsol04),
      tkIsolHEEP03CalcCfg_(tkIsolHEEP03),
      tkIsolHEEP04CalcCfg_(tkIsolHEEP04),
      magneticFieldToken_{cc.esConsumes()},
      caloGeometryToken_{cc.esConsumes()},
      caloTopologyToken_{cc.esConsumes()},
      trackerGeometryToken_{cc.esConsumes()},
      ecalSeveretyLevelAlgoToken_{cc.esConsumes()},
      ecalPFRechitThresholdsToken_{cc.esConsumes()},
      hcalHelperCone_{hcalCone, std::move(cc)},
      hcalHelperBc_{hcalBc, std::move(cc)},
      crackCorrectionFunction_{std::forward<std::unique_ptr<EcalClusterFunctionBaseClass>>(crackCorrectionFunction)},
      regHelper_{reg, cfg_.strategy.useEcalRegression, cfg_.strategy.useCombinationRegression, cc}

{
  ecalisoAlgo_ = std::make_unique<ElectronEcalPFClusterIsolation>(pfiso.ecaldrMax,
                                                                  pfiso.ecaldrVetoBarrel,
                                                                  pfiso.ecaldrVetoEndcap,
                                                                  pfiso.ecaletaStripBarrel,
                                                                  pfiso.ecaletaStripEndcap,
                                                                  pfiso.ecalenergyBarrel,
                                                                  pfiso.ecalenergyEndcap);
  hcalisoAlgo_ = std::make_unique<ElectronHcalPFClusterIsolation>(pfiso.hcaldrMax,
                                                                  pfiso.hcaldrVetoBarrel,
                                                                  pfiso.hcaldrVetoEndcap,
                                                                  pfiso.hcaletaStripBarrel,
                                                                  pfiso.hcaletaStripEndcap,
                                                                  pfiso.hcalenergyBarrel,
                                                                  pfiso.hcalenergyEndcap,
                                                                  pfiso.hcaluseEt);
}

void GsfElectronAlgo::checkSetup(const edm::EventSetup& es) {
  if (cfg_.strategy.useEcalRegression || cfg_.strategy.useCombinationRegression)
    regHelper_.checkSetup(es);

  if (crackCorrectionFunction_) {
    crackCorrectionFunction_->init(es);
  }
}

GsfElectronAlgo::EventData GsfElectronAlgo::beginEvent(edm::Event const& event,
                                                       CaloGeometry const& caloGeometry,
                                                       EcalSeverityLevelAlgo const& ecalSeveretyLevelAlgo) {
  auto const& hbheRecHits = event.get(cfg_.tokens.hbheRecHitsTag);

  // Isolation algos
  float egHcalIsoConeSizeOutSmall = 0.3, egHcalIsoConeSizeOutLarge = 0.4;
  float egHcalIsoConeSizeIn = cfg_.iso.intRadiusHcal, egHcalIsoPtMin = cfg_.iso.etMinHcal;

  float egIsoConeSizeOutSmall = 0.3, egIsoConeSizeOutLarge = 0.4, egIsoJurassicWidth = cfg_.iso.jurassicWidth;
  float egIsoPtMinBarrel = cfg_.iso.etMinBarrel, egIsoEMinBarrel = cfg_.iso.eMinBarrel,
        egIsoConeSizeInBarrel = cfg_.iso.intRadiusEcalBarrel;
  float egIsoPtMinEndcap = cfg_.iso.etMinEndcaps, egIsoEMinEndcap = cfg_.iso.eMinEndcaps,
        egIsoConeSizeInEndcap = cfg_.iso.intRadiusEcalEndcaps;

  auto barrelRecHits = event.getHandle(cfg_.tokens.barrelRecHitCollection);
  auto endcapRecHits = event.getHandle(cfg_.tokens.endcapRecHitCollection);

  edm::Handle<reco::PFClusterCollection> ecalPFClusters;
  std::vector<edm::Handle<reco::PFClusterCollection>> hcalPFClusters;
  if (cfg_.strategy.computePfClusterIso) {
    ecalPFClusters = event.getHandle(cfg_.tokens.pfClusterProducer);
    hcalPFClusters.push_back(event.getHandle(cfg_.tokens.pfClusterProducerHCAL));
    if (cfg_.pfiso.useHF) {
      hcalPFClusters.push_back(event.getHandle(cfg_.tokens.pfClusterProducerHFEM));
      hcalPFClusters.push_back(event.getHandle(cfg_.tokens.pfClusterProducerHFHAD));
    }
  }

  auto ctfTracks = event.getHandle(cfg_.tokens.ctfTracks);

  EventData eventData{
      .event = &event,
      .beamspot = &event.get(cfg_.tokens.beamSpotTag),
      .coreElectrons = event.getHandle(cfg_.tokens.gsfElectronCores),
      .barrelRecHits = barrelRecHits,
      .endcapRecHits = endcapRecHits,
      .currentCtfTracks = ctfTracks,
      .seeds = event.getHandle(cfg_.tokens.seedsTag),
      .vertices = event.getHandle(cfg_.tokens.vtxCollectionTag),
      .conversions = cfg_.strategy.fillConvVtxFitProb ? event.getHandle(cfg_.tokens.conversions)
                                                      : edm::Handle<reco::ConversionCollection>(),

      .hadIsolation03 = EgammaHcalIsolation(
          EgammaHcalIsolation::InclusionRule::withinConeAroundCluster,
          egHcalIsoConeSizeOutSmall,
          EgammaHcalIsolation::InclusionRule::withinConeAroundCluster,
          egHcalIsoConeSizeIn,
          hcalHelperCone_.eThresHB(),
          EgammaHcalIsolation::arrayHB{{egHcalIsoPtMin, egHcalIsoPtMin, egHcalIsoPtMin, egHcalIsoPtMin}},
          hcalHelperCone_.maxSeverityHB(),
          hcalHelperCone_.eThresHE(),
          EgammaHcalIsolation::arrayHE{{egHcalIsoPtMin,
                                        egHcalIsoPtMin,
                                        egHcalIsoPtMin,
                                        egHcalIsoPtMin,
                                        egHcalIsoPtMin,
                                        egHcalIsoPtMin,
                                        egHcalIsoPtMin}},
          hcalHelperCone_.maxSeverityHE(),
          hbheRecHits,
          caloGeometry,
          *hcalHelperCone_.hcalTopology(),
          *hcalHelperCone_.hcalChannelQuality(),
          *hcalHelperCone_.hcalSevLvlComputer(),
          *hcalHelperCone_.towerMap()),
      .hadIsolation04 = EgammaHcalIsolation(
          EgammaHcalIsolation::InclusionRule::withinConeAroundCluster,
          egHcalIsoConeSizeOutLarge,
          EgammaHcalIsolation::InclusionRule::withinConeAroundCluster,
          egHcalIsoConeSizeIn,
          hcalHelperCone_.eThresHB(),
          EgammaHcalIsolation::arrayHB{{egHcalIsoPtMin, egHcalIsoPtMin, egHcalIsoPtMin, egHcalIsoPtMin}},
          hcalHelperCone_.maxSeverityHB(),
          hcalHelperCone_.eThresHE(),
          EgammaHcalIsolation::arrayHE{{egHcalIsoPtMin,
                                        egHcalIsoPtMin,
                                        egHcalIsoPtMin,
                                        egHcalIsoPtMin,
                                        egHcalIsoPtMin,
                                        egHcalIsoPtMin,
                                        egHcalIsoPtMin}},
          hcalHelperCone_.maxSeverityHE(),
          hbheRecHits,
          caloGeometry,
          *hcalHelperCone_.hcalTopology(),
          *hcalHelperCone_.hcalChannelQuality(),
          *hcalHelperCone_.hcalSevLvlComputer(),
          *hcalHelperCone_.towerMap()),
      .hadIsolation03Bc = EgammaHcalIsolation(
          EgammaHcalIsolation::InclusionRule::withinConeAroundCluster,
          egHcalIsoConeSizeOutSmall,
          EgammaHcalIsolation::InclusionRule::isBehindClusterSeed,
          0.,
          hcalHelperCone_.eThresHB(),
          EgammaHcalIsolation::arrayHB{{egHcalIsoPtMin, egHcalIsoPtMin, egHcalIsoPtMin, egHcalIsoPtMin}},
          hcalHelperCone_.maxSeverityHB(),
          hcalHelperCone_.eThresHE(),
          EgammaHcalIsolation::arrayHE{{egHcalIsoPtMin,
                                        egHcalIsoPtMin,
                                        egHcalIsoPtMin,
                                        egHcalIsoPtMin,
                                        egHcalIsoPtMin,
                                        egHcalIsoPtMin,
                                        egHcalIsoPtMin}},
          hcalHelperCone_.maxSeverityHE(),
          hbheRecHits,
          caloGeometry,
          *hcalHelperCone_.hcalTopology(),
          *hcalHelperCone_.hcalChannelQuality(),
          *hcalHelperCone_.hcalSevLvlComputer(),
          *hcalHelperCone_.towerMap()),
      .hadIsolation04Bc = EgammaHcalIsolation(
          EgammaHcalIsolation::InclusionRule::withinConeAroundCluster,
          egHcalIsoConeSizeOutLarge,
          EgammaHcalIsolation::InclusionRule::isBehindClusterSeed,
          0.,
          hcalHelperCone_.eThresHB(),
          EgammaHcalIsolation::arrayHB{{egHcalIsoPtMin, egHcalIsoPtMin, egHcalIsoPtMin, egHcalIsoPtMin}},
          hcalHelperCone_.maxSeverityHB(),
          hcalHelperCone_.eThresHE(),
          EgammaHcalIsolation::arrayHE{{egHcalIsoPtMin,
                                        egHcalIsoPtMin,
                                        egHcalIsoPtMin,
                                        egHcalIsoPtMin,
                                        egHcalIsoPtMin,
                                        egHcalIsoPtMin,
                                        egHcalIsoPtMin}},
          hcalHelperCone_.maxSeverityHE(),
          hbheRecHits,
          caloGeometry,
          *hcalHelperCone_.hcalTopology(),
          *hcalHelperCone_.hcalChannelQuality(),
          *hcalHelperCone_.hcalSevLvlComputer(),
          *hcalHelperCone_.towerMap()),

      .ecalBarrelIsol03 = EgammaRecHitIsolation(egIsoConeSizeOutSmall,
                                                egIsoConeSizeInBarrel,
                                                egIsoJurassicWidth,
                                                egIsoPtMinBarrel,
                                                egIsoEMinBarrel,
                                                &caloGeometry,
                                                *barrelRecHits,
                                                &ecalSeveretyLevelAlgo,
                                                DetId::Ecal),
      .ecalBarrelIsol04 = EgammaRecHitIsolation(egIsoConeSizeOutLarge,
                                                egIsoConeSizeInBarrel,
                                                egIsoJurassicWidth,
                                                egIsoPtMinBarrel,
                                                egIsoEMinBarrel,
                                                &caloGeometry,
                                                *barrelRecHits,
                                                &ecalSeveretyLevelAlgo,
                                                DetId::Ecal),
      .ecalEndcapIsol03 = EgammaRecHitIsolation(egIsoConeSizeOutSmall,
                                                egIsoConeSizeInEndcap,
                                                egIsoJurassicWidth,
                                                egIsoPtMinEndcap,
                                                egIsoEMinEndcap,
                                                &caloGeometry,
                                                *endcapRecHits,
                                                &ecalSeveretyLevelAlgo,
                                                DetId::Ecal),
      .ecalEndcapIsol04 = EgammaRecHitIsolation(egIsoConeSizeOutLarge,
                                                egIsoConeSizeInEndcap,
                                                egIsoJurassicWidth,
                                                egIsoPtMinEndcap,
                                                egIsoEMinEndcap,
                                                &caloGeometry,
                                                *endcapRecHits,
                                                &ecalSeveretyLevelAlgo,
                                                DetId::Ecal),
      .tkIsol03Calc = EleTkIsolFromCands(tkIsol03CalcCfg_, *ctfTracks),
      .tkIsol04Calc = EleTkIsolFromCands(tkIsol04CalcCfg_, *ctfTracks),
      .tkIsolHEEP03Calc = EleTkIsolFromCands(tkIsolHEEP03CalcCfg_, *ctfTracks),
      .tkIsolHEEP04Calc = EleTkIsolFromCands(tkIsolHEEP04CalcCfg_, *ctfTracks),
      .originalCtfTracks = {},
      .originalGsfTracks = {},

      .ecalClustersHandle = ecalPFClusters,
      .hcalClustersHandle = hcalPFClusters

  };

  eventData.ecalBarrelIsol03.setUseNumCrystals(cfg_.iso.useNumCrystals);
  eventData.ecalBarrelIsol03.setVetoClustered(cfg_.iso.vetoClustered);
  eventData.ecalBarrelIsol03.doSeverityChecks(eventData.barrelRecHits.product(),
                                              cfg_.recHits.recHitSeverityToBeExcludedBarrel);
  eventData.ecalBarrelIsol03.doFlagChecks(cfg_.recHits.recHitFlagsToBeExcludedBarrel);
  eventData.ecalBarrelIsol04.setUseNumCrystals(cfg_.iso.useNumCrystals);
  eventData.ecalBarrelIsol04.setVetoClustered(cfg_.iso.vetoClustered);
  eventData.ecalBarrelIsol04.doSeverityChecks(eventData.barrelRecHits.product(),
                                              cfg_.recHits.recHitSeverityToBeExcludedBarrel);
  eventData.ecalBarrelIsol04.doFlagChecks(cfg_.recHits.recHitFlagsToBeExcludedBarrel);
  eventData.ecalEndcapIsol03.setUseNumCrystals(cfg_.iso.useNumCrystals);
  eventData.ecalEndcapIsol03.setVetoClustered(cfg_.iso.vetoClustered);
  eventData.ecalEndcapIsol03.doSeverityChecks(eventData.endcapRecHits.product(),
                                              cfg_.recHits.recHitSeverityToBeExcludedEndcaps);
  eventData.ecalEndcapIsol03.doFlagChecks(cfg_.recHits.recHitFlagsToBeExcludedEndcaps);
  eventData.ecalEndcapIsol04.setUseNumCrystals(cfg_.iso.useNumCrystals);
  eventData.ecalEndcapIsol04.setVetoClustered(cfg_.iso.vetoClustered);
  eventData.ecalEndcapIsol04.doSeverityChecks(eventData.endcapRecHits.product(),
                                              cfg_.recHits.recHitSeverityToBeExcludedEndcaps);
  eventData.ecalEndcapIsol04.doFlagChecks(cfg_.recHits.recHitFlagsToBeExcludedEndcaps);

  return eventData;
}

reco::GsfElectronCollection GsfElectronAlgo::completeElectrons(edm::Event const& event,
                                                               edm::EventSetup const& eventSetup,
                                                               const GsfElectronAlgo::HeavyObjectCache* hoc,
                                                               const HcalPFCuts* hcalCuts) {
  reco::GsfElectronCollection electrons;

  auto const& magneticField = eventSetup.getData(magneticFieldToken_);
  auto const& caloGeometry = eventSetup.getData(caloGeometryToken_);
  auto const& caloTopology = eventSetup.getData(caloTopologyToken_);
  auto const& trackerGeometry = eventSetup.getData(trackerGeometryToken_);
  auto const& ecalSeveretyLevelAlgo = eventSetup.getData(ecalSeveretyLevelAlgoToken_);
  auto const& thresholds = eventSetup.getData(ecalPFRechitThresholdsToken_);

  // prepare access to hcal data
  hcalHelperCone_.beginEvent(event, eventSetup);
  hcalHelperBc_.beginEvent(event, eventSetup);

  checkSetup(eventSetup);
  auto eventData = beginEvent(event, caloGeometry, ecalSeveretyLevelAlgo);
  double magneticFieldInTesla = magneticField.inTesla(GlobalPoint(0., 0., 0.)).z();

  MultiTrajectoryStateTransform mtsTransform(&trackerGeometry, &magneticField);
  GsfConstraintAtVertex constraintAtVtx(&trackerGeometry, &magneticField);

  std::optional<egamma::conv::TrackTable> ctfTrackTable = std::nullopt;
  std::optional<egamma::conv::TrackTable> gsfTrackTable = std::nullopt;

  const GsfElectronCoreCollection* coreCollection = eventData.coreElectrons.product();
  for (unsigned int i = 0; i < coreCollection->size(); ++i) {
    // check there is no existing electron with this core
    const GsfElectronCoreRef coreRef = edm::Ref<GsfElectronCoreCollection>(eventData.coreElectrons, i);

    // check there is a super-cluster
    if (coreRef->superCluster().isNull())
      continue;

    // prepare internal structure for electron specific data
    ElectronData electronData(coreRef, *eventData.beamspot);

    // calculate and check Trajectory StatesOnSurface....
    if (!electronData.calculateTSOS(mtsTransform, constraintAtVtx))
      continue;

    eventData.retreiveOriginalTrackCollections(electronData.ctfTrackRef, electronData.coreRef->gsfTrack());

    if (!eventData.originalCtfTracks.isValid()) {
      eventData.originalCtfTracks = eventData.currentCtfTracks;
    }

    if (ctfTrackTable == std::nullopt) {
      ctfTrackTable = egamma::conv::TrackTable(*eventData.originalCtfTracks);
    }
    if (gsfTrackTable == std::nullopt) {
      gsfTrackTable = egamma::conv::TrackTable(*eventData.originalGsfTracks);
    }

    createElectron(electrons,
                   electronData,
                   eventData,
                   caloTopology,
                   caloGeometry,
                   mtsTransform,
                   magneticFieldInTesla,
                   hoc,
                   ctfTrackTable.value(),
                   gsfTrackTable.value(),
                   thresholds,
                   hcalCuts);

  }  // loop over tracks
  return electrons;
}

void GsfElectronAlgo::setCutBasedPreselectionFlag(GsfElectron& ele, const reco::BeamSpot& bs) const {
  // default value
  ele.setPassCutBasedPreselection(false);

  // kind of seeding
  bool eg = ele.core()->ecalDrivenSeed();
  bool pf = ele.core()->trackerDrivenSeed() && !ele.core()->ecalDrivenSeed();
  if (eg && pf) {
    throw cms::Exception("GsfElectronAlgo|BothEcalAndPureTrackerDriven")
        << "An electron cannot be both egamma and purely pflow";
  }
  if ((!eg) && (!pf)) {
    throw cms::Exception("GsfElectronAlgo|NeitherEcalNorPureTrackerDriven")
        << "An electron cannot be neither egamma nor purely pflow";
  }

  CutsConfiguration const& cfg = cfg_.cuts;

  // Et cut
  double etaValue = EleRelPoint(ele.superCluster()->position(), bs.position()).eta();
  double etValue = ele.superCluster()->energy() / cosh(etaValue);
  LogTrace("GsfElectronAlgo") << "Et : " << etValue;
  if (ele.isEB() && (etValue < cfg.minSCEtBarrel))
    return;
  if (ele.isEE() && (etValue < cfg.minSCEtEndcaps))
    return;
  LogTrace("GsfElectronAlgo") << "Et criteria are satisfied";

  // E/p cut
  double eopValue = ele.eSuperClusterOverP();
  LogTrace("GsfElectronAlgo") << "E/p : " << eopValue;
  if (ele.isEB() && (eopValue > cfg.maxEOverPBarrel))
    return;
  if (ele.isEE() && (eopValue > cfg.maxEOverPEndcaps))
    return;
  if (ele.isEB() && (eopValue < cfg.minEOverPBarrel))
    return;
  if (ele.isEE() && (eopValue < cfg.minEOverPEndcaps))
    return;
  LogTrace("GsfElectronAlgo") << "E/p criteria are satisfied";

  // HoE cuts
  LogTrace("GsfElectronAlgo") << "HoE : " << ele.hcalOverEcal();
  double hoeCone = ele.hcalOverEcal();
  double hoeBc = ele.hcalOverEcalBc();
  const reco::CaloCluster& seedCluster = *(ele.superCluster()->seed());
  int detector = seedCluster.hitsAndFractions()[0].first.subdetId();
  bool HoEveto = false;
  double scle = ele.superCluster()->energy();

  if (detector == EcalBarrel)
    HoEveto = hoeCone * scle < cfg.maxHBarrelCone || hoeBc * scle < cfg.maxHBarrelBc ||
              hoeCone < cfg.maxHOverEBarrelCone || hoeBc < cfg.maxHOverEBarrelBc;
  else if (detector == EcalEndcap)
    HoEveto = hoeCone * scle < cfg.maxHEndcapsCone || hoeBc * scle < cfg.maxHEndcapsBc ||
              hoeCone < cfg.maxHOverEEndcapsCone || hoeBc < cfg.maxHOverEEndcapsBc;

  if (!HoEveto)
    return;
  LogTrace("GsfElectronAlgo") << "H/E criteria are satisfied";

  // delta eta criteria
  double deta = ele.deltaEtaSuperClusterTrackAtVtx();
  LogTrace("GsfElectronAlgo") << "delta eta : " << deta;
  if (ele.isEB() && (std::abs(deta) > cfg.maxDeltaEtaBarrel))
    return;
  if (ele.isEE() && (std::abs(deta) > cfg.maxDeltaEtaEndcaps))
    return;
  LogTrace("GsfElectronAlgo") << "Delta eta criteria are satisfied";

  // delta phi criteria
  double dphi = ele.deltaPhiSuperClusterTrackAtVtx();
  LogTrace("GsfElectronAlgo") << "delta phi : " << dphi;
  if (ele.isEB() && (std::abs(dphi) > cfg.maxDeltaPhiBarrel))
    return;
  if (ele.isEE() && (std::abs(dphi) > cfg.maxDeltaPhiEndcaps))
    return;
  LogTrace("GsfElectronAlgo") << "Delta phi criteria are satisfied";

  // sigma ieta ieta
  LogTrace("GsfElectronAlgo") << "sigma ieta ieta : " << ele.sigmaIetaIeta();
  if (ele.isEB() && (ele.sigmaIetaIeta() > cfg.maxSigmaIetaIetaBarrel))
    return;
  if (ele.isEE() && (ele.sigmaIetaIeta() > cfg.maxSigmaIetaIetaEndcaps))
    return;
  LogTrace("GsfElectronAlgo") << "Sigma ieta ieta criteria are satisfied";

  // fiducial
  if (!ele.isEB() && cfg.isBarrel)
    return;
  if (!ele.isEE() && cfg.isEndcaps)
    return;
  if (cfg.isFiducial &&
      (ele.isEBEEGap() || ele.isEBEtaGap() || ele.isEBPhiGap() || ele.isEERingGap() || ele.isEEDeeGap()))
    return;
  LogTrace("GsfElectronAlgo") << "Fiducial flags criteria are satisfied";

  // seed in TEC
  edm::RefToBase<TrajectorySeed> seed = ele.gsfTrack()->extra()->seedRef();
  ElectronSeedRef elseed = seed.castTo<ElectronSeedRef>();
  if (eg && !cfg_.cuts.seedFromTEC) {
    if (elseed.isNull()) {
      throw cms::Exception("GsfElectronAlgo|NotElectronSeed") << "The GsfTrack seed is not an ElectronSeed ?!";
    } else {
      if (elseed->subDet(1) == 6)
        return;
    }
  }

  // transverse impact parameter
  if (std::abs(ele.gsfTrack()->dxy(bs.position())) > cfg.maxTIP)
    return;
  LogTrace("GsfElectronAlgo") << "TIP criterion is satisfied";

  LogTrace("GsfElectronAlgo") << "All cut based criteria are satisfied";
  ele.setPassCutBasedPreselection(true);
}

void GsfElectronAlgo::createElectron(reco::GsfElectronCollection& electrons,
                                     ElectronData& electronData,
                                     EventData& eventData,
                                     CaloTopology const& topology,
                                     CaloGeometry const& geometry,
                                     MultiTrajectoryStateTransform const& mtsTransform,
                                     double magneticFieldInTesla,
                                     const GsfElectronAlgo::HeavyObjectCache* hoc,
                                     egamma::conv::TrackTableView ctfTable,
                                     egamma::conv::TrackTableView gsfTable,
                                     EcalPFRecHitThresholds const& thresholds,
                                     const HcalPFCuts* hcalCuts) {
  // charge ID
  int eleCharge;
  GsfElectron::ChargeInfo eleChargeInfo;
  electronData.computeCharge(eleCharge, eleChargeInfo);

  // electron basic cluster
  CaloClusterPtr elbcRef = electronData.getEleBasicCluster(mtsTransform);

  // Seed cluster
  const reco::CaloCluster& seedCluster = *(electronData.superClusterRef->seed());

  // seed Xtal
  // temporary, till CaloCluster->seed() is made available
  DetId seedXtalId = seedCluster.hitsAndFractions()[0].first;

  electronData.calculateMode();

  //====================================================
  // Candidate attributes
  //====================================================

  Candidate::LorentzVector momentum = electronData.calculateMomentum();

  //====================================================
  // Track-Cluster Matching
  //====================================================

  reco::GsfElectron::TrackClusterMatching tcMatching;
  tcMatching.electronCluster = elbcRef;
  tcMatching.eSuperClusterOverP =
      (electronData.vtxMom.mag() > 0) ? (electronData.superClusterRef->energy() / electronData.vtxMom.mag()) : (-1.);
  tcMatching.eSeedClusterOverP =
      (electronData.vtxMom.mag() > 0.) ? (seedCluster.energy() / electronData.vtxMom.mag()) : (-1);
  tcMatching.eSeedClusterOverPout =
      (electronData.seedMom.mag() > 0.) ? (seedCluster.energy() / electronData.seedMom.mag()) : (-1.);
  tcMatching.eEleClusterOverPout =
      (electronData.eleMom.mag() > 0.) ? (elbcRef->energy() / electronData.eleMom.mag()) : (-1.);

  EleRelPointPair scAtVtx(
      electronData.superClusterRef->position(), electronData.sclPos, eventData.beamspot->position());
  tcMatching.deltaEtaSuperClusterAtVtx = scAtVtx.dEta();
  tcMatching.deltaPhiSuperClusterAtVtx = scAtVtx.dPhi();

  EleRelPointPair seedAtCalo(seedCluster.position(), electronData.seedPos, eventData.beamspot->position());
  tcMatching.deltaEtaSeedClusterAtCalo = seedAtCalo.dEta();
  tcMatching.deltaPhiSeedClusterAtCalo = seedAtCalo.dPhi();

  EleRelPointPair ecAtCalo(elbcRef->position(), electronData.elePos, eventData.beamspot->position());
  tcMatching.deltaEtaEleClusterAtCalo = ecAtCalo.dEta();
  tcMatching.deltaPhiEleClusterAtCalo = ecAtCalo.dPhi();

  //=======================================================
  // Track extrapolations
  //=======================================================

  reco::GsfElectron::TrackExtrapolations tkExtra;
  ele_convert(electronData.vtxPos, tkExtra.positionAtVtx);
  ele_convert(electronData.sclPos, tkExtra.positionAtCalo);
  ele_convert(electronData.vtxMom, tkExtra.momentumAtVtx);
  ele_convert(electronData.sclMom, tkExtra.momentumAtCalo);
  ele_convert(electronData.seedMom, tkExtra.momentumOut);
  ele_convert(electronData.eleMom, tkExtra.momentumAtEleClus);
  ele_convert(electronData.vtxMomWithConstraint, tkExtra.momentumAtVtxWithConstraint);

  //=======================================================
  // Closest Ctf Track
  //=======================================================

  reco::GsfElectron::ClosestCtfTrack ctfInfo;
  ctfInfo.ctfTrack = electronData.ctfTrackRef;
  ctfInfo.shFracInnerHits = electronData.shFracInnerHits;

  //====================================================
  // FiducialFlags, using nextToBoundary definition of gaps
  //====================================================

  reco::GsfElectron::FiducialFlags fiducialFlags;
  int region = seedXtalId.det();
  int detector = seedXtalId.subdetId();
  double feta = std::abs(electronData.superClusterRef->position().eta());
  if (detector == EcalBarrel) {
    fiducialFlags.isEB = true;
    EBDetId ebdetid(seedXtalId);
    if (EBDetId::isNextToEtaBoundary(ebdetid)) {
      if (ebdetid.ietaAbs() == 85) {
        fiducialFlags.isEBEEGap = true;
      } else {
        fiducialFlags.isEBEtaGap = true;
      }
    }
    if (EBDetId::isNextToPhiBoundary(ebdetid)) {
      fiducialFlags.isEBPhiGap = true;
    }
  } else if (detector == EcalEndcap) {
    fiducialFlags.isEE = true;
    EEDetId eedetid(seedXtalId);
    if (EEDetId::isNextToRingBoundary(eedetid)) {
      if (std::abs(feta) < 2.) {
        fiducialFlags.isEBEEGap = true;
      } else {
        fiducialFlags.isEERingGap = true;
      }
    }
    if (EEDetId::isNextToDBoundary(eedetid)) {
      fiducialFlags.isEEDeeGap = true;
    }
  } else if (EcalTools::isHGCalDet((DetId::Detector)region)) {
    fiducialFlags.isEE = true;
    //HGCalDetId eeDetid(seedXtalId);
    // fill in fiducial information when we know how to use it...
  } else {
    throw cms::Exception("GsfElectronAlgo|UnknownXtalRegion")
        << "createElectron(): do not know if it is a barrel or endcap seed cluster !!!!";
  }

  //====================================================
  // SaturationInfo
  //====================================================

  auto saturationInfo = calculateSaturationInfo(electronData.superClusterRef, eventData);

  //====================================================
  // ShowerShape
  //====================================================

  reco::GsfElectron::ShowerShape showerShape;
  reco::GsfElectron::ShowerShape full5x5_showerShape;
  if (!EcalTools::isHGCalDet((DetId::Detector)region)) {
    showerShape = calculateShowerShape<false>(electronData.superClusterRef,
                                              hcalHelperCone_,
                                              hcalHelperBc_,
                                              eventData,
                                              topology,
                                              geometry,
                                              thresholds,
                                              hcalCuts);
    full5x5_showerShape = calculateShowerShape<true>(electronData.superClusterRef,
                                                     hcalHelperCone_,
                                                     hcalHelperBc_,
                                                     eventData,
                                                     topology,
                                                     geometry,
                                                     thresholds,
                                                     hcalCuts);
  }

  //====================================================
  // ConversionRejection
  //====================================================

  edm::Handle<reco::TrackCollection> ctfTracks = eventData.originalCtfTracks;
  if (!ctfTracks.isValid()) {
    ctfTracks = eventData.currentCtfTracks;
  }

  {
    //get the references to the gsf and ctf tracks that are made
    //by the electron
    const reco::TrackRef el_ctftrack = electronData.coreRef->ctfTrack();
    const reco::GsfTrackRef& el_gsftrack = electronData.coreRef->gsfTrack();

    //protect against the wrong collection being passed to the function
    if (el_ctftrack.isNonnull() && el_ctftrack.id() != ctfTracks.id())
      throw cms::Exception("ConversionFinderError")
          << "ProductID of ctf track collection does not match ProductID of electron's CTF track! \n";
    if (el_gsftrack.isNonnull() && el_gsftrack.id() != eventData.originalGsfTracks.id())
      throw cms::Exception("ConversionFinderError")
          << "ProductID of gsf track collection does not match ProductID of electron's GSF track! \n";
  }

  // values of conversionInfo.flag
  // -9999 : Partner track was not found
  // 0     : Partner track found in the CTF collection using
  // 1     : Partner track found in the CTF collection using
  // 2     : Partner track found in the GSF collection using
  // 3     : Partner track found in the GSF collection using the electron's GSF track
  auto conversionInfo = egamma::conv::findConversion(*electronData.coreRef, ctfTable, gsfTable, magneticFieldInTesla);

  reco::GsfElectron::ConversionRejection conversionVars;
  conversionVars.flags = conversionInfo.flag;
  conversionVars.dist = conversionInfo.dist;
  conversionVars.dcot = conversionInfo.dcot;
  conversionVars.radius = conversionInfo.radiusOfConversion;
  if (cfg_.strategy.fillConvVtxFitProb) {
    //this is an intentionally bugged version which ignores the GsfTrack
    //this is a bug which was introduced in reduced e/gamma where the GsfTrack gets
    //relinked to a new collection which means it can no longer match the conversion
    //as it matches based on product/id
    //we keep this defination for the MVAs
    const auto matchedConv = ConversionTools::matchedConversion(
        electronData.coreRef->ctfTrack(), *eventData.conversions, eventData.beamspot->position(), 2.0, 1e-6, 0);
    conversionVars.vtxFitProb = ConversionTools::getVtxFitProb(matchedConv);
  }
  if (conversionInfo.conversionPartnerCtfTkIdx) {
    conversionVars.partner = TrackBaseRef(reco::TrackRef(ctfTracks, conversionInfo.conversionPartnerCtfTkIdx.value()));
  } else if (conversionInfo.conversionPartnerGsfTkIdx) {
    conversionVars.partner =
        TrackBaseRef(reco::GsfTrackRef(eventData.originalGsfTracks, conversionInfo.conversionPartnerGsfTkIdx.value()));
  }

  //====================================================
  // Go !
  //====================================================

  electrons.emplace_back(eleCharge,
                         eleChargeInfo,
                         electronData.coreRef,
                         tcMatching,
                         tkExtra,
                         ctfInfo,
                         fiducialFlags,
                         showerShape,
                         full5x5_showerShape,
                         conversionVars,
                         saturationInfo);
  auto& ele = electrons.back();
  // Will be overwritten later in the case of the regression
  ele.setCorrectedEcalEnergyError(egamma::ecalClusterEnergyUncertaintyElectronSpecific(*(ele.superCluster())));
  ele.setP4(GsfElectron::P4_FROM_SUPER_CLUSTER, momentum, 0, true);
  ele.setMass(0.0);

  //====================================================
  // brems fractions
  //====================================================

  if (electronData.innMom.mag() > 0.) {
    ele.setTrackFbrem((electronData.innMom.mag() - electronData.outMom.mag()) / electronData.innMom.mag());
  }

  // the supercluster is the refined one The seed is not necessarily the first cluster
  // hence the use of the electronCluster
  SuperClusterRef sc = ele.superCluster();
  if (!(sc.isNull())) {
    CaloClusterPtr cl = ele.electronCluster();
    if (sc->clustersSize() > 1) {
      float pf_fbrem = (sc->energy() - cl->energy()) / sc->energy();
      ele.setSuperClusterFbrem(pf_fbrem);
    } else {
      ele.setSuperClusterFbrem(0);
    }
  }

  //====================================================
  // classification and corrections
  //====================================================
  // classification
  const auto elClass = egamma::classifyElectron(ele);
  ele.setClassification(elClass);

  bool unexpectedClassification = elClass == GsfElectron::UNKNOWN || elClass > GsfElectron::GAP;
  if (unexpectedClassification) {
    edm::LogWarning("GsfElectronAlgo") << "unexpected classification";
  }

  // ecal energy
  if (cfg_.strategy.useEcalRegression)  // new
  {
    regHelper_.applyEcalRegression(ele, *eventData.vertices, *eventData.barrelRecHits, *eventData.endcapRecHits);
  } else  // original implementation
  {
    if (!EcalTools::isHGCalDet((DetId::Detector)region)) {
      if (ele.core()->ecalDrivenSeed()) {
        if (cfg_.strategy.ecalDrivenEcalEnergyFromClassBasedParameterization && !unexpectedClassification) {
          if (ele.isEcalEnergyCorrected()) {
            edm::LogWarning("ElectronEnergyCorrector::classBasedElectronEnergy") << "already done";
          } else {
            ele.setCorrectedEcalEnergy(
                egamma::classBasedElectronEnergy(ele, *eventData.beamspot, *crackCorrectionFunction_));
          }
        }
        if (cfg_.strategy.ecalDrivenEcalErrorFromClassBasedParameterization) {
          ele.setCorrectedEcalEnergyError(egamma::classBasedElectronEnergyUncertainty(ele));
        }
      } else {
        if (cfg_.strategy.pureTrackerDrivenEcalErrorFromSimpleParameterization) {
          ele.setCorrectedEcalEnergyError(egamma::simpleElectronEnergyUncertainty(ele));
        }
      }
    }
  }

  // momentum
  // Keep the default correction running first. The track momentum error is computed in there
  if (cfg_.strategy.useDefaultEnergyCorrection && ele.core()->ecalDrivenSeed() && !unexpectedClassification) {
    if (ele.p4Error(reco::GsfElectron::P4_COMBINATION) != 999.) {
      edm::LogWarning("ElectronMomentumCorrector::correct") << "already done";
    } else {
      auto p = egamma::correctElectronMomentum(ele, electronData.vtxTSOS);
      ele.correctMomentum(p.momentum, p.trackError, p.finalError);
    }
  }
  if (cfg_.strategy.useCombinationRegression)  // new
  {
    regHelper_.applyCombinationRegression(ele);
  }

  //====================================================
  // now isolation variables
  //====================================================
  reco::GsfElectron::IsolationVariables dr03, dr04;
  dr03.tkSumPt = eventData.tkIsol03Calc(*ele.gsfTrack()).ptSum;
  dr04.tkSumPt = eventData.tkIsol04Calc(*ele.gsfTrack()).ptSum;
  dr03.tkSumPtHEEP = eventData.tkIsolHEEP03Calc(*ele.gsfTrack()).ptSum;
  dr04.tkSumPtHEEP = eventData.tkIsolHEEP04Calc(*ele.gsfTrack()).ptSum;

  if (!EcalTools::isHGCalDet((DetId::Detector)region)) {
    for (uint id = 0; id < dr03.hcalRecHitSumEt.size(); ++id) {
      dr03.hcalRecHitSumEt[id] = eventData.hadIsolation03.getHcalEtSum(&ele, id + 1, hcalCuts);
      dr03.hcalRecHitSumEtBc[id] = eventData.hadIsolation03Bc.getHcalEtSumBc(&ele, id + 1, hcalCuts);

      dr04.hcalRecHitSumEt[id] = eventData.hadIsolation04.getHcalEtSum(&ele, id + 1, hcalCuts);
      dr04.hcalRecHitSumEtBc[id] = eventData.hadIsolation04Bc.getHcalEtSumBc(&ele, id + 1, hcalCuts);
    }

    dr03.ecalRecHitSumEt = eventData.ecalBarrelIsol03.getEtSum(&ele, thresholds);
    dr03.ecalRecHitSumEt += eventData.ecalEndcapIsol03.getEtSum(&ele, thresholds);

    dr04.ecalRecHitSumEt = eventData.ecalBarrelIsol04.getEtSum(&ele, thresholds);
    dr04.ecalRecHitSumEt += eventData.ecalEndcapIsol04.getEtSum(&ele, thresholds);
  }

  dr03.pre7DepthHcal = false;
  dr04.pre7DepthHcal = false;

  ele.setIsolation03(dr03);
  ele.setIsolation04(dr04);

  //====================================================
  // PFclusters based ISO !
  // Added in CMSSW_12_0_1 at this stage to be able to use them in PF ID DNN
  //====================================================
  if (cfg_.strategy.computePfClusterIso) {
    reco::GsfElectron::PflowIsolationVariables isoVariables;
    isoVariables.sumEcalClusterEt = ecalisoAlgo_->getSum(ele, eventData.ecalClustersHandle);
    isoVariables.sumHcalClusterEt = hcalisoAlgo_->getSum(ele, eventData.hcalClustersHandle);
    // Other Pfiso variables are initialized at 0 and not used
    ele.setPfIsolationVariables(isoVariables);
  }

  //====================================================
  // preselection flag
  //====================================================

  setCutBasedPreselectionFlag(ele, *eventData.beamspot);
  //setting mva flag, currently GedGsfElectron and GsfElectron pre-selection flags have desynced
  //this is for GedGsfElectrons, GsfElectrons (ie old pre 7X std reco) resets this later on
  //in the function "addPfInfo"
  //yes this is awful, we'll fix it once we work out how to...
  float mvaValue = hoc->sElectronMVAEstimator->mva(ele, *(eventData.vertices));
  ele.setPassMvaPreselection(mvaValue > cfg_.strategy.PreSelectMVA);

  //====================================================
  // Pixel match variables
  //====================================================
  setPixelMatchInfomation(ele);

  LogTrace("GsfElectronAlgo") << "Constructed new electron with energy  " << ele.p4().e();
}

// Pixel match variables
void GsfElectronAlgo::setPixelMatchInfomation(reco::GsfElectron& ele) const {
  int sd1 = 0;
  int sd2 = 0;
  float dPhi1 = 0;
  float dPhi2 = 0;
  float dRz1 = 0;
  float dRz2 = 0;
  edm::RefToBase<TrajectorySeed> seed = ele.gsfTrack()->extra()->seedRef();
  ElectronSeedRef elseed = seed.castTo<ElectronSeedRef>();
  if (seed.isNull()) {
  } else {
    if (elseed.isNull()) {
    } else {
      sd1 = elseed->subDet(0);
      sd2 = elseed->subDet(1);
      dPhi1 = (ele.charge() > 0) ? elseed->dPhiPos(0) : elseed->dPhiNeg(0);
      dPhi2 = (ele.charge() > 0) ? elseed->dPhiPos(1) : elseed->dPhiNeg(1);
      dRz1 = (ele.charge() > 0) ? elseed->dRZPos(0) : elseed->dRZNeg(0);
      dRz2 = (ele.charge() > 0) ? elseed->dRZPos(1) : elseed->dRZNeg(1);
    }
  }
  ele.setPixelMatchSubdetectors(sd1, sd2);
  ele.setPixelMatchDPhi1(dPhi1);
  ele.setPixelMatchDPhi2(dPhi2);
  ele.setPixelMatchDRz1(dRz1);
  ele.setPixelMatchDRz2(dRz2);
}