Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:24:48 

0001 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0002 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0003 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0004 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
0005 #include "DataFormats/EgammaCandidates/interface/GsfElectronFwd.h"
0006 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
0007 #include "DataFormats/EgammaReco/interface/ElectronSeed.h"
0008 #include "DataFormats/ForwardDetId/interface/ForwardSubdetector.h"
0009 #include "DataFormats/ForwardDetId/interface/HGCalDetId.h"
0010 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0011 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0012 #include "DataFormats/Math/interface/LorentzVector.h"
0013 #include "DataFormats/ParticleFlowReco/interface/GsfPFRecTrack.h"
0014 #include "DataFormats/TrackingRecHit/interface/TrackingRecHitFwd.h"
0015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0016 #include "FWCore/Utilities/interface/isFinite.h"
0017 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0018 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h"
0019 #include "RecoEcal/EgammaCoreTools/interface/EcalTools.h"
0020 #include "RecoEgamma/EgammaElectronAlgos/interface/EgAmbiguityTools.h"
0021 #include "RecoEgamma/EgammaElectronAlgos/interface/ElectronClassification.h"
0022 #include "RecoEgamma/EgammaElectronAlgos/interface/ElectronEnergyCorrector.h"
0023 #include "RecoEgamma/EgammaElectronAlgos/interface/ElectronMomentumCorrector.h"
0024 #include "RecoEgamma/EgammaElectronAlgos/interface/ElectronUtilities.h"
0025 #include "RecoEgamma/EgammaElectronAlgos/interface/GsfElectronAlgo.h"
0026 #include "RecoEgamma/EgammaElectronAlgos/interface/ecalClusterEnergyUncertaintyElectronSpecific.h"
0027 #include "CommonTools/Egamma/interface/ConversionTools.h"
0028 #include "RecoEcal/EgammaCoreTools/interface/EgammaLocalCovParamDefaults.h"
0029 #include "CondFormats/DataRecord/interface/HcalPFCutsRcd.h"
0030 #include "CondTools/Hcal/interface/HcalPFCutsHandler.h"
0031 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0032 
0033 #include <Math/Point3D.h>
0034 #include <memory>
0035 
0036 #include <algorithm>
0037 #include <sstream>
0038 
0039 using namespace edm;
0040 using namespace reco;
0041 
0042 GsfElectronAlgo::HeavyObjectCache::HeavyObjectCache(const edm::ParameterSet& conf) {
0043   // soft electron MVA
0044   SoftElectronMVAEstimator::Configuration sconfig;
0045   sconfig.vweightsfiles = conf.getParameter<std::vector<std::string>>("SoftElecMVAFilesString");
0046   sElectronMVAEstimator = std::make_unique<SoftElectronMVAEstimator>(sconfig);
0047   // isolated electron MVA
0048   ElectronMVAEstimator::Configuration iconfig;
0049   iconfig.vweightsfiles = conf.getParameter<std::vector<std::string>>("ElecMVAFilesString");
0050   iElectronMVAEstimator = std::make_unique<ElectronMVAEstimator>(iconfig);
0051 
0052   // Here we will have to load the DNN PFID if present in the config
0053   egammaTools::DNNConfiguration dconfig;
0054   const auto& pset_dnn = conf.getParameter<edm::ParameterSet>("EleDNNPFid");
0055   const bool dnnEnabled = pset_dnn.getParameter<bool>("enabled");
0056   if (dnnEnabled) {
0057     dconfig.inputTensorName = pset_dnn.getParameter<std::string>("inputTensorName");
0058     dconfig.outputTensorName = pset_dnn.getParameter<std::string>("outputTensorName");
0059     dconfig.modelsFiles = pset_dnn.getParameter<std::vector<std::string>>("modelsFiles");
0060     dconfig.scalersFiles = pset_dnn.getParameter<std::vector<std::string>>("scalersFiles");
0061     dconfig.outputDim = pset_dnn.getParameter<std::vector<unsigned int>>("outputDim");
0062     const auto useEBModelInGap = pset_dnn.getParameter<bool>("useEBModelInGap");
0063     iElectronDNNEstimator = std::make_unique<ElectronDNNEstimator>(dconfig, useEBModelInGap);
0064   }
0065 }
0066 
0067 //===================================================================
0068 // GsfElectronAlgo::EventData
0069 //===================================================================
0070 
0071 struct GsfElectronAlgo::EventData {
0072   // utilities
0073   void retreiveOriginalTrackCollections(const reco::TrackRef&, const reco::GsfTrackRef&);
0074 
0075   // general
0076   edm::Event const* event;
0077   const reco::BeamSpot* beamspot;
0078 
0079   // input collections
0080   edm::Handle<reco::GsfElectronCoreCollection> coreElectrons;
0081   edm::Handle<EcalRecHitCollection> barrelRecHits;
0082   edm::Handle<EcalRecHitCollection> endcapRecHits;
0083   edm::Handle<reco::TrackCollection> currentCtfTracks;
0084   edm::Handle<reco::ElectronSeedCollection> seeds;
0085   edm::Handle<reco::VertexCollection> vertices;
0086   edm::Handle<reco::ConversionCollection> conversions;
0087 
0088   // isolation helpers
0089   EgammaHcalIsolation hadIsolation03, hadIsolation04;
0090   EgammaHcalIsolation hadIsolation03Bc, hadIsolation04Bc;
0091   EgammaRecHitIsolation ecalBarrelIsol03, ecalBarrelIsol04;
0092   EgammaRecHitIsolation ecalEndcapIsol03, ecalEndcapIsol04;
0093 
0094   EleTkIsolFromCands tkIsol03Calc;
0095   EleTkIsolFromCands tkIsol04Calc;
0096   EleTkIsolFromCands tkIsolHEEP03Calc;
0097   EleTkIsolFromCands tkIsolHEEP04Calc;
0098 
0099   edm::Handle<reco::TrackCollection> originalCtfTracks;
0100   edm::Handle<reco::GsfTrackCollection> originalGsfTracks;
0101 
0102   bool originalCtfTrackCollectionRetreived = false;
0103   bool originalGsfTrackCollectionRetreived = false;
0104 
0105   //PFCluster Isolation handles
0106   edm::Handle<reco::PFClusterCollection> ecalClustersHandle;
0107   std::vector<edm::Handle<reco::PFClusterCollection>> hcalClustersHandle;
0108 };
0109 
0110 //===================================================================
0111 // GsfElectronAlgo::ElectronData
0112 //===================================================================
0113 
0114 struct GsfElectronAlgo::ElectronData {
0115   // Refs to subproducts
0116   const reco::GsfElectronCoreRef coreRef;
0117   const reco::GsfTrackRef gsfTrackRef;
0118   const reco::SuperClusterRef superClusterRef;
0119   reco::TrackRef ctfTrackRef;
0120   float shFracInnerHits;
0121   const reco::BeamSpot beamSpot;
0122 
0123   // constructors
0124   ElectronData(const reco::GsfElectronCoreRef& core, const reco::BeamSpot& bs);
0125 
0126   // utilities
0127   void computeCharge(int& charge, reco::GsfElectron::ChargeInfo& info);
0128   reco::CaloClusterPtr getEleBasicCluster(MultiTrajectoryStateTransform const&);
0129   bool calculateTSOS(MultiTrajectoryStateTransform const&, GsfConstraintAtVertex const&);
0130   void calculateMode();
0131   reco::Candidate::LorentzVector calculateMomentum();
0132 
0133   // TSOS
0134   TrajectoryStateOnSurface innTSOS;
0135   TrajectoryStateOnSurface outTSOS;
0136   TrajectoryStateOnSurface vtxTSOS;
0137   TrajectoryStateOnSurface sclTSOS;
0138   TrajectoryStateOnSurface seedTSOS;
0139   TrajectoryStateOnSurface eleTSOS;
0140   TrajectoryStateOnSurface constrainedVtxTSOS;
0141 
0142   // mode
0143   GlobalVector innMom, seedMom, eleMom, sclMom, vtxMom, outMom;
0144   GlobalPoint innPos, seedPos, elePos, sclPos, vtxPos, outPos;
0145   GlobalVector vtxMomWithConstraint;
0146 };
0147 
0148 void GsfElectronAlgo::EventData::retreiveOriginalTrackCollections(const reco::TrackRef& ctfTrack,
0149                                                                   const reco::GsfTrackRef& gsfTrack) {
0150   if ((!originalCtfTrackCollectionRetreived) && (ctfTrack.isNonnull())) {
0151     event->get(ctfTrack.id(), originalCtfTracks);
0152     originalCtfTrackCollectionRetreived = true;
0153   }
0154   if ((!originalGsfTrackCollectionRetreived) && (gsfTrack.isNonnull())) {
0155     event->get(gsfTrack.id(), originalGsfTracks);
0156     originalGsfTrackCollectionRetreived = true;
0157   }
0158 }
0159 
0160 GsfElectronAlgo::ElectronData::ElectronData(const reco::GsfElectronCoreRef& core, const reco::BeamSpot& bs)
0161     : coreRef(core),
0162       gsfTrackRef(coreRef->gsfTrack()),
0163       superClusterRef(coreRef->superCluster()),
0164       ctfTrackRef(coreRef->ctfTrack()),
0165       shFracInnerHits(coreRef->ctfGsfOverlap()),
0166       beamSpot(bs) {}
0167 
0168 void GsfElectronAlgo::ElectronData::computeCharge(int& charge, GsfElectron::ChargeInfo& info) {
0169   // determine charge from SC
0170   GlobalPoint orig, scpos;
0171   ele_convert(beamSpot.position(), orig);
0172   ele_convert(superClusterRef->position(), scpos);
0173   GlobalVector scvect(scpos - orig);
0174   GlobalPoint inntkpos = innTSOS.globalPosition();
0175   GlobalVector inntkvect = GlobalVector(inntkpos - orig);
0176   float dPhiInnEle = normalizedPhi(scvect.barePhi() - inntkvect.barePhi());
0177   if (dPhiInnEle > 0)
0178     info.scPixCharge = -1;
0179   else
0180     info.scPixCharge = 1;
0181 
0182   // flags
0183   int chargeGsf = gsfTrackRef->charge();
0184   info.isGsfScPixConsistent = ((chargeGsf * info.scPixCharge) > 0);
0185   info.isGsfCtfConsistent = (ctfTrackRef.isNonnull() && ((chargeGsf * ctfTrackRef->charge()) > 0));
0186   info.isGsfCtfScPixConsistent = (info.isGsfScPixConsistent && info.isGsfCtfConsistent);
0187 
0188   // default charge
0189   if (info.isGsfScPixConsistent || ctfTrackRef.isNull()) {
0190     charge = info.scPixCharge;
0191   } else {
0192     charge = ctfTrackRef->charge();
0193   }
0194 }
0195 
0196 CaloClusterPtr GsfElectronAlgo::ElectronData::getEleBasicCluster(MultiTrajectoryStateTransform const& mtsTransform) {
0197   CaloClusterPtr eleRef;
0198   TrajectoryStateOnSurface tempTSOS;
0199   TrajectoryStateOnSurface outTSOS = mtsTransform.outerStateOnSurface(*gsfTrackRef);
0200   float dphimin = 1.e30;
0201   for (auto const& bc : superClusterRef->clusters()) {
0202     GlobalPoint posclu(bc->position().x(), bc->position().y(), bc->position().z());
0203     tempTSOS = mtsTransform.extrapolatedState(outTSOS, posclu);
0204     if (!tempTSOS.isValid())
0205       tempTSOS = outTSOS;
0206     GlobalPoint extrap = tempTSOS.globalPosition();
0207     float dphi = EleRelPointPair(posclu, extrap, beamSpot.position()).dPhi();
0208     if (std::abs(dphi) < dphimin) {
0209       dphimin = std::abs(dphi);
0210       eleRef = bc;
0211       eleTSOS = tempTSOS;
0212     }
0213   }
0214   return eleRef;
0215 }
0216 
0217 bool GsfElectronAlgo::ElectronData::calculateTSOS(MultiTrajectoryStateTransform const& mtsTransform,
0218                                                   GsfConstraintAtVertex const& constraintAtVtx) {
0219   //at innermost point
0220   innTSOS = mtsTransform.innerStateOnSurface(*gsfTrackRef);
0221   if (!innTSOS.isValid())
0222     return false;
0223 
0224   //at vertex
0225   // innermost state propagation to the beam spot position
0226   GlobalPoint bsPos;
0227   ele_convert(beamSpot.position(), bsPos);
0228   vtxTSOS = mtsTransform.extrapolatedState(innTSOS, bsPos);
0229   if (!vtxTSOS.isValid())
0230     vtxTSOS = innTSOS;
0231 
0232   //at seed
0233   outTSOS = mtsTransform.outerStateOnSurface(*gsfTrackRef);
0234   if (!outTSOS.isValid())
0235     return false;
0236 
0237   //    TrajectoryStateOnSurface seedTSOS
0238   seedTSOS = mtsTransform.extrapolatedState(outTSOS,
0239                                             GlobalPoint(superClusterRef->seed()->position().x(),
0240                                                         superClusterRef->seed()->position().y(),
0241                                                         superClusterRef->seed()->position().z()));
0242   if (!seedTSOS.isValid())
0243     seedTSOS = outTSOS;
0244 
0245   // at scl
0246   sclTSOS = mtsTransform.extrapolatedState(
0247       innTSOS, GlobalPoint(superClusterRef->x(), superClusterRef->y(), superClusterRef->z()));
0248   if (!sclTSOS.isValid())
0249     sclTSOS = outTSOS;
0250 
0251   // constrained momentum
0252   constrainedVtxTSOS = constraintAtVtx.constrainAtBeamSpot(*gsfTrackRef, beamSpot);
0253 
0254   return true;
0255 }
0256 
0257 void GsfElectronAlgo::ElectronData::calculateMode() {
0258   multiTrajectoryStateMode::momentumFromModeCartesian(innTSOS, innMom);
0259   multiTrajectoryStateMode::positionFromModeCartesian(innTSOS, innPos);
0260   multiTrajectoryStateMode::momentumFromModeCartesian(seedTSOS, seedMom);
0261   multiTrajectoryStateMode::positionFromModeCartesian(seedTSOS, seedPos);
0262   multiTrajectoryStateMode::momentumFromModeCartesian(eleTSOS, eleMom);
0263   multiTrajectoryStateMode::positionFromModeCartesian(eleTSOS, elePos);
0264   multiTrajectoryStateMode::momentumFromModeCartesian(sclTSOS, sclMom);
0265   multiTrajectoryStateMode::positionFromModeCartesian(sclTSOS, sclPos);
0266   multiTrajectoryStateMode::momentumFromModeCartesian(vtxTSOS, vtxMom);
0267   multiTrajectoryStateMode::positionFromModeCartesian(vtxTSOS, vtxPos);
0268   multiTrajectoryStateMode::momentumFromModeCartesian(outTSOS, outMom);
0269   multiTrajectoryStateMode::positionFromModeCartesian(outTSOS, outPos);
0270   multiTrajectoryStateMode::momentumFromModeCartesian(constrainedVtxTSOS, vtxMomWithConstraint);
0271 }
0272 
0273 Candidate::LorentzVector GsfElectronAlgo::ElectronData::calculateMomentum() {
0274   double scale = superClusterRef->energy() / vtxMom.mag();
0275   return Candidate::LorentzVector(
0276       vtxMom.x() * scale, vtxMom.y() * scale, vtxMom.z() * scale, superClusterRef->energy());
0277 }
0278 
0279 reco::GsfElectron::SaturationInfo GsfElectronAlgo::calculateSaturationInfo(const reco::SuperClusterRef& theClus,
0280                                                                            EventData const& eventData) const {
0281   reco::GsfElectron::SaturationInfo si;
0282 
0283   const reco::CaloCluster& seedCluster = *(theClus->seed());
0284   DetId seedXtalId = seedCluster.seed();
0285   int detector = seedXtalId.subdetId();
0286   const EcalRecHitCollection* ecalRecHits = nullptr;
0287   if (detector == EcalBarrel)
0288     ecalRecHits = eventData.barrelRecHits.product();
0289   else
0290     ecalRecHits = eventData.endcapRecHits.product();
0291 
0292   int nSaturatedXtals = 0;
0293   bool isSeedSaturated = false;
0294   for (auto&& hitFractionPair : theClus->hitsAndFractions()) {
0295     auto&& ecalRecHit = ecalRecHits->find(hitFractionPair.first);
0296     if (ecalRecHit == ecalRecHits->end())
0297       continue;
0298     if (ecalRecHit->checkFlag(EcalRecHit::Flags::kSaturated)) {
0299       nSaturatedXtals++;
0300       if (seedXtalId == ecalRecHit->detid())
0301         isSeedSaturated = true;
0302     }
0303   }
0304   si.nSaturatedXtals = nSaturatedXtals;
0305   si.isSeedSaturated = isSeedSaturated;
0306 
0307   return si;
0308 }
0309 
0310 template <bool full5x5>
0311 reco::GsfElectron::ShowerShape GsfElectronAlgo::calculateShowerShape(const reco::SuperClusterRef& theClus,
0312                                                                      ElectronHcalHelper const& hcalHelperCone,
0313                                                                      ElectronHcalHelper const& hcalHelperBc,
0314                                                                      EventData const& eventData,
0315                                                                      CaloTopology const& topology,
0316                                                                      CaloGeometry const& geometry,
0317                                                                      EcalPFRecHitThresholds const& thresholds,
0318                                                                      const HcalPFCuts* hcalCuts) const {
0319   using ClusterTools = EcalClusterToolsT<full5x5>;
0320   reco::GsfElectron::ShowerShape showerShape;
0321 
0322   const reco::CaloCluster& seedCluster = *(theClus->seed());
0323   // temporary, till CaloCluster->seed() is made available
0324   DetId seedXtalId = seedCluster.hitsAndFractions()[0].first;
0325   int detector = seedXtalId.subdetId();
0326 
0327   const EcalRecHitCollection* recHits = nullptr;
0328   std::vector<int> recHitFlagsToBeExcluded;
0329   std::vector<int> recHitSeverityToBeExcluded;
0330   if (detector == EcalBarrel) {
0331     recHits = eventData.barrelRecHits.product();
0332     recHitFlagsToBeExcluded = cfg_.recHits.recHitFlagsToBeExcludedBarrel;
0333     recHitSeverityToBeExcluded = cfg_.recHits.recHitSeverityToBeExcludedBarrel;
0334   } else {
0335     recHits = eventData.endcapRecHits.product();
0336     recHitFlagsToBeExcluded = cfg_.recHits.recHitFlagsToBeExcludedEndcaps;
0337     recHitSeverityToBeExcluded = cfg_.recHits.recHitSeverityToBeExcludedEndcaps;
0338   }
0339 
0340   const auto& covariances = ClusterTools::covariances(seedCluster, recHits, &topology, &geometry);
0341 
0342   // do noise-cleaning for full5x5, by passing per crystal PF recHit thresholds and mult values
0343   // mult values for EB and EE were obtained by dedicated studies
0344   const auto& localCovariances = full5x5 ? ClusterTools::localCovariances(seedCluster,
0345                                                                           recHits,
0346                                                                           &topology,
0347                                                                           EgammaLocalCovParamDefaults::kRelEnCut,
0348                                                                           &thresholds,
0349                                                                           cfg_.cuts.multThresEB,
0350                                                                           cfg_.cuts.multThresEE)
0351                                          : ClusterTools::localCovariances(seedCluster, recHits, &topology);
0352 
0353   showerShape.sigmaEtaEta = sqrt(covariances[0]);
0354   showerShape.sigmaIetaIeta = sqrt(localCovariances[0]);
0355   if (!edm::isNotFinite(localCovariances[2]))
0356     showerShape.sigmaIphiIphi = sqrt(localCovariances[2]);
0357   showerShape.e1x5 = ClusterTools::e1x5(seedCluster, recHits, &topology);
0358   showerShape.e2x5Max = ClusterTools::e2x5Max(seedCluster, recHits, &topology);
0359   showerShape.e5x5 = ClusterTools::e5x5(seedCluster, recHits, &topology);
0360   showerShape.r9 = ClusterTools::e3x3(seedCluster, recHits, &topology) / theClus->rawEnergy();
0361 
0362   const float scale = full5x5 ? showerShape.e5x5 : theClus->energy();
0363 
0364   for (uint id = 0; id < showerShape.hcalOverEcal.size(); ++id) {
0365     showerShape.hcalOverEcal[id] = hcalHelperCone.hcalESum(*theClus, id + 1, hcalCuts) / scale;
0366     showerShape.hcalOverEcalBc[id] = hcalHelperBc.hcalESum(*theClus, id + 1, hcalCuts) / scale;
0367   }
0368   showerShape.invalidHcal = !hcalHelperBc.hasActiveHcal(*theClus);
0369   showerShape.hcalTowersBehindClusters = hcalHelperBc.hcalTowersBehindClusters(*theClus);
0370   showerShape.pre7DepthHcal = false;
0371 
0372   // extra shower shapes
0373   const float see_by_spp = showerShape.sigmaIetaIeta * showerShape.sigmaIphiIphi;
0374   if (see_by_spp > 0) {
0375     showerShape.sigmaIetaIphi = localCovariances[1] / see_by_spp;
0376   } else if (localCovariances[1] > 0) {
0377     showerShape.sigmaIetaIphi = 1.f;
0378   } else {
0379     showerShape.sigmaIetaIphi = -1.f;
0380   }
0381   showerShape.eMax = ClusterTools::eMax(seedCluster, recHits);
0382   showerShape.e2nd = ClusterTools::e2nd(seedCluster, recHits);
0383   showerShape.eTop = ClusterTools::eTop(seedCluster, recHits, &topology);
0384   showerShape.eLeft = ClusterTools::eLeft(seedCluster, recHits, &topology);
0385   showerShape.eRight = ClusterTools::eRight(seedCluster, recHits, &topology);
0386   showerShape.eBottom = ClusterTools::eBottom(seedCluster, recHits, &topology);
0387 
0388   showerShape.e2x5Left = ClusterTools::e2x5Left(seedCluster, recHits, &topology);
0389   showerShape.e2x5Right = ClusterTools::e2x5Right(seedCluster, recHits, &topology);
0390   showerShape.e2x5Top = ClusterTools::e2x5Top(seedCluster, recHits, &topology);
0391   showerShape.e2x5Bottom = ClusterTools::e2x5Bottom(seedCluster, recHits, &topology);
0392 
0393   return showerShape;
0394 }
0395 
0396 //===================================================================
0397 // GsfElectronAlgo
0398 //===================================================================
0399 
0400 GsfElectronAlgo::GsfElectronAlgo(const Tokens& input,
0401                                  const StrategyConfiguration& strategy,
0402                                  const CutsConfiguration& cuts,
0403                                  const ElectronHcalHelper::Configuration& hcalCone,
0404                                  const ElectronHcalHelper::Configuration& hcalBc,
0405                                  const IsolationConfiguration& iso,
0406                                  const PFClusterIsolationConfiguration& pfiso,
0407                                  const EcalRecHitsConfiguration& recHits,
0408                                  std::unique_ptr<EcalClusterFunctionBaseClass>&& crackCorrectionFunction,
0409                                  const RegressionHelper::Configuration& reg,
0410                                  const edm::ParameterSet& tkIsol03,
0411                                  const edm::ParameterSet& tkIsol04,
0412                                  const edm::ParameterSet& tkIsolHEEP03,
0413                                  const edm::ParameterSet& tkIsolHEEP04,
0414                                  edm::ConsumesCollector&& cc)
0415     : cfg_{input, strategy, cuts, iso, pfiso, recHits},
0416       tkIsol03CalcCfg_(tkIsol03),
0417       tkIsol04CalcCfg_(tkIsol04),
0418       tkIsolHEEP03CalcCfg_(tkIsolHEEP03),
0419       tkIsolHEEP04CalcCfg_(tkIsolHEEP04),
0420       magneticFieldToken_{cc.esConsumes()},
0421       caloGeometryToken_{cc.esConsumes()},
0422       caloTopologyToken_{cc.esConsumes()},
0423       trackerGeometryToken_{cc.esConsumes()},
0424       ecalSeveretyLevelAlgoToken_{cc.esConsumes()},
0425       ecalPFRechitThresholdsToken_{cc.esConsumes()},
0426       hcalHelperCone_{hcalCone, std::move(cc)},
0427       hcalHelperBc_{hcalBc, std::move(cc)},
0428       crackCorrectionFunction_{std::forward<std::unique_ptr<EcalClusterFunctionBaseClass>>(crackCorrectionFunction)},
0429       regHelper_{reg, cfg_.strategy.useEcalRegression, cfg_.strategy.useCombinationRegression, cc}
0430 
0431 {
0432   ///PF ECAL cluster based isolations
0433   ecalisoAlgo_ = std::make_unique<ElectronEcalPFClusterIsolation>(pfiso.ecaldrMax,
0434                                                                   pfiso.ecaldrVetoBarrel,
0435                                                                   pfiso.ecaldrVetoEndcap,
0436                                                                   pfiso.ecaletaStripBarrel,
0437                                                                   pfiso.ecaletaStripEndcap,
0438                                                                   pfiso.ecalenergyBarrel,
0439                                                                   pfiso.ecalenergyEndcap);
0440   hcalisoAlgo_ = std::make_unique<ElectronHcalPFClusterIsolation>(pfiso.hcaldrMax,
0441                                                                   pfiso.hcaldrVetoBarrel,
0442                                                                   pfiso.hcaldrVetoEndcap,
0443                                                                   pfiso.hcaletaStripBarrel,
0444                                                                   pfiso.hcaletaStripEndcap,
0445                                                                   pfiso.hcalenergyBarrel,
0446                                                                   pfiso.hcalenergyEndcap,
0447                                                                   pfiso.hcaluseEt);
0448 }
0449 
0450 void GsfElectronAlgo::checkSetup(const edm::EventSetup& es) {
0451   if (cfg_.strategy.useEcalRegression || cfg_.strategy.useCombinationRegression)
0452     regHelper_.checkSetup(es);
0453 
0454   if (crackCorrectionFunction_) {
0455     crackCorrectionFunction_->init(es);
0456   }
0457 }
0458 
0459 GsfElectronAlgo::EventData GsfElectronAlgo::beginEvent(edm::Event const& event,
0460                                                        CaloGeometry const& caloGeometry,
0461                                                        EcalSeverityLevelAlgo const& ecalSeveretyLevelAlgo) {
0462   auto const& hbheRecHits = event.get(cfg_.tokens.hbheRecHitsTag);
0463 
0464   // Isolation algos
0465   float egHcalIsoConeSizeOutSmall = 0.3, egHcalIsoConeSizeOutLarge = 0.4;
0466   float egHcalIsoConeSizeIn = cfg_.iso.intRadiusHcal, egHcalIsoPtMin = cfg_.iso.etMinHcal;
0467 
0468   float egIsoConeSizeOutSmall = 0.3, egIsoConeSizeOutLarge = 0.4, egIsoJurassicWidth = cfg_.iso.jurassicWidth;
0469   float egIsoPtMinBarrel = cfg_.iso.etMinBarrel, egIsoEMinBarrel = cfg_.iso.eMinBarrel,
0470         egIsoConeSizeInBarrel = cfg_.iso.intRadiusEcalBarrel;
0471   float egIsoPtMinEndcap = cfg_.iso.etMinEndcaps, egIsoEMinEndcap = cfg_.iso.eMinEndcaps,
0472         egIsoConeSizeInEndcap = cfg_.iso.intRadiusEcalEndcaps;
0473 
0474   auto barrelRecHits = event.getHandle(cfg_.tokens.barrelRecHitCollection);
0475   auto endcapRecHits = event.getHandle(cfg_.tokens.endcapRecHitCollection);
0476 
0477   edm::Handle<reco::PFClusterCollection> ecalPFClusters;
0478   std::vector<edm::Handle<reco::PFClusterCollection>> hcalPFClusters;
0479   if (cfg_.strategy.computePfClusterIso) {
0480     ecalPFClusters = event.getHandle(cfg_.tokens.pfClusterProducer);
0481     hcalPFClusters.push_back(event.getHandle(cfg_.tokens.pfClusterProducerHCAL));
0482     if (cfg_.pfiso.useHF) {
0483       hcalPFClusters.push_back(event.getHandle(cfg_.tokens.pfClusterProducerHFEM));
0484       hcalPFClusters.push_back(event.getHandle(cfg_.tokens.pfClusterProducerHFHAD));
0485     }
0486   }
0487 
0488   auto ctfTracks = event.getHandle(cfg_.tokens.ctfTracks);
0489 
0490   EventData eventData{
0491       .event = &event,
0492       .beamspot = &event.get(cfg_.tokens.beamSpotTag),
0493       .coreElectrons = event.getHandle(cfg_.tokens.gsfElectronCores),
0494       .barrelRecHits = barrelRecHits,
0495       .endcapRecHits = endcapRecHits,
0496       .currentCtfTracks = ctfTracks,
0497       .seeds = event.getHandle(cfg_.tokens.seedsTag),
0498       .vertices = event.getHandle(cfg_.tokens.vtxCollectionTag),
0499       .conversions = cfg_.strategy.fillConvVtxFitProb ? event.getHandle(cfg_.tokens.conversions)
0500                                                       : edm::Handle<reco::ConversionCollection>(),
0501 
0502       .hadIsolation03 = EgammaHcalIsolation(
0503           EgammaHcalIsolation::InclusionRule::withinConeAroundCluster,
0504           egHcalIsoConeSizeOutSmall,
0505           EgammaHcalIsolation::InclusionRule::withinConeAroundCluster,
0506           egHcalIsoConeSizeIn,
0507           hcalHelperCone_.eThresHB(),
0508           EgammaHcalIsolation::arrayHB{{egHcalIsoPtMin, egHcalIsoPtMin, egHcalIsoPtMin, egHcalIsoPtMin}},
0509           hcalHelperCone_.maxSeverityHB(),
0510           hcalHelperCone_.eThresHE(),
0511           EgammaHcalIsolation::arrayHE{{egHcalIsoPtMin,
0512                                         egHcalIsoPtMin,
0513                                         egHcalIsoPtMin,
0514                                         egHcalIsoPtMin,
0515                                         egHcalIsoPtMin,
0516                                         egHcalIsoPtMin,
0517                                         egHcalIsoPtMin}},
0518           hcalHelperCone_.maxSeverityHE(),
0519           hbheRecHits,
0520           caloGeometry,
0521           *hcalHelperCone_.hcalTopology(),
0522           *hcalHelperCone_.hcalChannelQuality(),
0523           *hcalHelperCone_.hcalSevLvlComputer(),
0524           *hcalHelperCone_.towerMap()),
0525       .hadIsolation04 = EgammaHcalIsolation(
0526           EgammaHcalIsolation::InclusionRule::withinConeAroundCluster,
0527           egHcalIsoConeSizeOutLarge,
0528           EgammaHcalIsolation::InclusionRule::withinConeAroundCluster,
0529           egHcalIsoConeSizeIn,
0530           hcalHelperCone_.eThresHB(),
0531           EgammaHcalIsolation::arrayHB{{egHcalIsoPtMin, egHcalIsoPtMin, egHcalIsoPtMin, egHcalIsoPtMin}},
0532           hcalHelperCone_.maxSeverityHB(),
0533           hcalHelperCone_.eThresHE(),
0534           EgammaHcalIsolation::arrayHE{{egHcalIsoPtMin,
0535                                         egHcalIsoPtMin,
0536                                         egHcalIsoPtMin,
0537                                         egHcalIsoPtMin,
0538                                         egHcalIsoPtMin,
0539                                         egHcalIsoPtMin,
0540                                         egHcalIsoPtMin}},
0541           hcalHelperCone_.maxSeverityHE(),
0542           hbheRecHits,
0543           caloGeometry,
0544           *hcalHelperCone_.hcalTopology(),
0545           *hcalHelperCone_.hcalChannelQuality(),
0546           *hcalHelperCone_.hcalSevLvlComputer(),
0547           *hcalHelperCone_.towerMap()),
0548       .hadIsolation03Bc = EgammaHcalIsolation(
0549           EgammaHcalIsolation::InclusionRule::withinConeAroundCluster,
0550           egHcalIsoConeSizeOutSmall,
0551           EgammaHcalIsolation::InclusionRule::isBehindClusterSeed,
0552           0.,
0553           hcalHelperCone_.eThresHB(),
0554           EgammaHcalIsolation::arrayHB{{egHcalIsoPtMin, egHcalIsoPtMin, egHcalIsoPtMin, egHcalIsoPtMin}},
0555           hcalHelperCone_.maxSeverityHB(),
0556           hcalHelperCone_.eThresHE(),
0557           EgammaHcalIsolation::arrayHE{{egHcalIsoPtMin,
0558                                         egHcalIsoPtMin,
0559                                         egHcalIsoPtMin,
0560                                         egHcalIsoPtMin,
0561                                         egHcalIsoPtMin,
0562                                         egHcalIsoPtMin,
0563                                         egHcalIsoPtMin}},
0564           hcalHelperCone_.maxSeverityHE(),
0565           hbheRecHits,
0566           caloGeometry,
0567           *hcalHelperCone_.hcalTopology(),
0568           *hcalHelperCone_.hcalChannelQuality(),
0569           *hcalHelperCone_.hcalSevLvlComputer(),
0570           *hcalHelperCone_.towerMap()),
0571       .hadIsolation04Bc = EgammaHcalIsolation(
0572           EgammaHcalIsolation::InclusionRule::withinConeAroundCluster,
0573           egHcalIsoConeSizeOutLarge,
0574           EgammaHcalIsolation::InclusionRule::isBehindClusterSeed,
0575           0.,
0576           hcalHelperCone_.eThresHB(),
0577           EgammaHcalIsolation::arrayHB{{egHcalIsoPtMin, egHcalIsoPtMin, egHcalIsoPtMin, egHcalIsoPtMin}},
0578           hcalHelperCone_.maxSeverityHB(),
0579           hcalHelperCone_.eThresHE(),
0580           EgammaHcalIsolation::arrayHE{{egHcalIsoPtMin,
0581                                         egHcalIsoPtMin,
0582                                         egHcalIsoPtMin,
0583                                         egHcalIsoPtMin,
0584                                         egHcalIsoPtMin,
0585                                         egHcalIsoPtMin,
0586                                         egHcalIsoPtMin}},
0587           hcalHelperCone_.maxSeverityHE(),
0588           hbheRecHits,
0589           caloGeometry,
0590           *hcalHelperCone_.hcalTopology(),
0591           *hcalHelperCone_.hcalChannelQuality(),
0592           *hcalHelperCone_.hcalSevLvlComputer(),
0593           *hcalHelperCone_.towerMap()),
0594 
0595       .ecalBarrelIsol03 = EgammaRecHitIsolation(egIsoConeSizeOutSmall,
0596                                                 egIsoConeSizeInBarrel,
0597                                                 egIsoJurassicWidth,
0598                                                 egIsoPtMinBarrel,
0599                                                 egIsoEMinBarrel,
0600                                                 &caloGeometry,
0601                                                 *barrelRecHits,
0602                                                 &ecalSeveretyLevelAlgo,
0603                                                 DetId::Ecal),
0604       .ecalBarrelIsol04 = EgammaRecHitIsolation(egIsoConeSizeOutLarge,
0605                                                 egIsoConeSizeInBarrel,
0606                                                 egIsoJurassicWidth,
0607                                                 egIsoPtMinBarrel,
0608                                                 egIsoEMinBarrel,
0609                                                 &caloGeometry,
0610                                                 *barrelRecHits,
0611                                                 &ecalSeveretyLevelAlgo,
0612                                                 DetId::Ecal),
0613       .ecalEndcapIsol03 = EgammaRecHitIsolation(egIsoConeSizeOutSmall,
0614                                                 egIsoConeSizeInEndcap,
0615                                                 egIsoJurassicWidth,
0616                                                 egIsoPtMinEndcap,
0617                                                 egIsoEMinEndcap,
0618                                                 &caloGeometry,
0619                                                 *endcapRecHits,
0620                                                 &ecalSeveretyLevelAlgo,
0621                                                 DetId::Ecal),
0622       .ecalEndcapIsol04 = EgammaRecHitIsolation(egIsoConeSizeOutLarge,
0623                                                 egIsoConeSizeInEndcap,
0624                                                 egIsoJurassicWidth,
0625                                                 egIsoPtMinEndcap,
0626                                                 egIsoEMinEndcap,
0627                                                 &caloGeometry,
0628                                                 *endcapRecHits,
0629                                                 &ecalSeveretyLevelAlgo,
0630                                                 DetId::Ecal),
0631       .tkIsol03Calc = EleTkIsolFromCands(tkIsol03CalcCfg_, *ctfTracks),
0632       .tkIsol04Calc = EleTkIsolFromCands(tkIsol04CalcCfg_, *ctfTracks),
0633       .tkIsolHEEP03Calc = EleTkIsolFromCands(tkIsolHEEP03CalcCfg_, *ctfTracks),
0634       .tkIsolHEEP04Calc = EleTkIsolFromCands(tkIsolHEEP04CalcCfg_, *ctfTracks),
0635       .originalCtfTracks = {},
0636       .originalGsfTracks = {},
0637 
0638       .ecalClustersHandle = ecalPFClusters,
0639       .hcalClustersHandle = hcalPFClusters
0640 
0641   };
0642 
0643   eventData.ecalBarrelIsol03.setUseNumCrystals(cfg_.iso.useNumCrystals);
0644   eventData.ecalBarrelIsol03.setVetoClustered(cfg_.iso.vetoClustered);
0645   eventData.ecalBarrelIsol03.doSeverityChecks(eventData.barrelRecHits.product(),
0646                                               cfg_.recHits.recHitSeverityToBeExcludedBarrel);
0647   eventData.ecalBarrelIsol03.doFlagChecks(cfg_.recHits.recHitFlagsToBeExcludedBarrel);
0648   eventData.ecalBarrelIsol04.setUseNumCrystals(cfg_.iso.useNumCrystals);
0649   eventData.ecalBarrelIsol04.setVetoClustered(cfg_.iso.vetoClustered);
0650   eventData.ecalBarrelIsol04.doSeverityChecks(eventData.barrelRecHits.product(),
0651                                               cfg_.recHits.recHitSeverityToBeExcludedBarrel);
0652   eventData.ecalBarrelIsol04.doFlagChecks(cfg_.recHits.recHitFlagsToBeExcludedBarrel);
0653   eventData.ecalEndcapIsol03.setUseNumCrystals(cfg_.iso.useNumCrystals);
0654   eventData.ecalEndcapIsol03.setVetoClustered(cfg_.iso.vetoClustered);
0655   eventData.ecalEndcapIsol03.doSeverityChecks(eventData.endcapRecHits.product(),
0656                                               cfg_.recHits.recHitSeverityToBeExcludedEndcaps);
0657   eventData.ecalEndcapIsol03.doFlagChecks(cfg_.recHits.recHitFlagsToBeExcludedEndcaps);
0658   eventData.ecalEndcapIsol04.setUseNumCrystals(cfg_.iso.useNumCrystals);
0659   eventData.ecalEndcapIsol04.setVetoClustered(cfg_.iso.vetoClustered);
0660   eventData.ecalEndcapIsol04.doSeverityChecks(eventData.endcapRecHits.product(),
0661                                               cfg_.recHits.recHitSeverityToBeExcludedEndcaps);
0662   eventData.ecalEndcapIsol04.doFlagChecks(cfg_.recHits.recHitFlagsToBeExcludedEndcaps);
0663 
0664   return eventData;
0665 }
0666 
0667 reco::GsfElectronCollection GsfElectronAlgo::completeElectrons(edm::Event const& event,
0668                                                                edm::EventSetup const& eventSetup,
0669                                                                const GsfElectronAlgo::HeavyObjectCache* hoc,
0670                                                                const HcalPFCuts* hcalCuts) {
0671   reco::GsfElectronCollection electrons;
0672 
0673   auto const& magneticField = eventSetup.getData(magneticFieldToken_);
0674   auto const& caloGeometry = eventSetup.getData(caloGeometryToken_);
0675   auto const& caloTopology = eventSetup.getData(caloTopologyToken_);
0676   auto const& trackerGeometry = eventSetup.getData(trackerGeometryToken_);
0677   auto const& ecalSeveretyLevelAlgo = eventSetup.getData(ecalSeveretyLevelAlgoToken_);
0678   auto const& thresholds = eventSetup.getData(ecalPFRechitThresholdsToken_);
0679 
0680   // prepare access to hcal data
0681   hcalHelperCone_.beginEvent(event, eventSetup);
0682   hcalHelperBc_.beginEvent(event, eventSetup);
0683 
0684   checkSetup(eventSetup);
0685   auto eventData = beginEvent(event, caloGeometry, ecalSeveretyLevelAlgo);
0686   double magneticFieldInTesla = magneticField.inTesla(GlobalPoint(0., 0., 0.)).z();
0687 
0688   MultiTrajectoryStateTransform mtsTransform(&trackerGeometry, &magneticField);
0689   GsfConstraintAtVertex constraintAtVtx(&trackerGeometry, &magneticField);
0690 
0691   std::optional<egamma::conv::TrackTable> ctfTrackTable = std::nullopt;
0692   std::optional<egamma::conv::TrackTable> gsfTrackTable = std::nullopt;
0693 
0694   const GsfElectronCoreCollection* coreCollection = eventData.coreElectrons.product();
0695   for (unsigned int i = 0; i < coreCollection->size(); ++i) {
0696     // check there is no existing electron with this core
0697     const GsfElectronCoreRef coreRef = edm::Ref<GsfElectronCoreCollection>(eventData.coreElectrons, i);
0698 
0699     // check there is a super-cluster
0700     if (coreRef->superCluster().isNull())
0701       continue;
0702 
0703     // prepare internal structure for electron specific data
0704     ElectronData electronData(coreRef, *eventData.beamspot);
0705 
0706     // calculate and check Trajectory StatesOnSurface....
0707     if (!electronData.calculateTSOS(mtsTransform, constraintAtVtx))
0708       continue;
0709 
0710     eventData.retreiveOriginalTrackCollections(electronData.ctfTrackRef, electronData.coreRef->gsfTrack());
0711 
0712     if (!eventData.originalCtfTracks.isValid()) {
0713       eventData.originalCtfTracks = eventData.currentCtfTracks;
0714     }
0715 
0716     if (ctfTrackTable == std::nullopt) {
0717       ctfTrackTable = egamma::conv::TrackTable(*eventData.originalCtfTracks);
0718     }
0719     if (gsfTrackTable == std::nullopt) {
0720       gsfTrackTable = egamma::conv::TrackTable(*eventData.originalGsfTracks);
0721     }
0722 
0723     createElectron(electrons,
0724                    electronData,
0725                    eventData,
0726                    caloTopology,
0727                    caloGeometry,
0728                    mtsTransform,
0729                    magneticFieldInTesla,
0730                    hoc,
0731                    ctfTrackTable.value(),
0732                    gsfTrackTable.value(),
0733                    thresholds,
0734                    hcalCuts);
0735 
0736   }  // loop over tracks
0737   return electrons;
0738 }
0739 
0740 void GsfElectronAlgo::setCutBasedPreselectionFlag(GsfElectron& ele, const reco::BeamSpot& bs) const {
0741   // default value
0742   ele.setPassCutBasedPreselection(false);
0743 
0744   // kind of seeding
0745   bool eg = ele.core()->ecalDrivenSeed();
0746   bool pf = ele.core()->trackerDrivenSeed() && !ele.core()->ecalDrivenSeed();
0747   if (eg && pf) {
0748     throw cms::Exception("GsfElectronAlgo|BothEcalAndPureTrackerDriven")
0749         << "An electron cannot be both egamma and purely pflow";
0750   }
0751   if ((!eg) && (!pf)) {
0752     throw cms::Exception("GsfElectronAlgo|NeitherEcalNorPureTrackerDriven")
0753         << "An electron cannot be neither egamma nor purely pflow";
0754   }
0755 
0756   CutsConfiguration const& cfg = cfg_.cuts;
0757 
0758   // Et cut
0759   double etaValue = EleRelPoint(ele.superCluster()->position(), bs.position()).eta();
0760   double etValue = ele.superCluster()->energy() / cosh(etaValue);
0761   LogTrace("GsfElectronAlgo") << "Et : " << etValue;
0762   if (ele.isEB() && (etValue < cfg.minSCEtBarrel))
0763     return;
0764   if (ele.isEE() && (etValue < cfg.minSCEtEndcaps))
0765     return;
0766   LogTrace("GsfElectronAlgo") << "Et criteria are satisfied";
0767 
0768   // E/p cut
0769   double eopValue = ele.eSuperClusterOverP();
0770   LogTrace("GsfElectronAlgo") << "E/p : " << eopValue;
0771   if (ele.isEB() && (eopValue > cfg.maxEOverPBarrel))
0772     return;
0773   if (ele.isEE() && (eopValue > cfg.maxEOverPEndcaps))
0774     return;
0775   if (ele.isEB() && (eopValue < cfg.minEOverPBarrel))
0776     return;
0777   if (ele.isEE() && (eopValue < cfg.minEOverPEndcaps))
0778     return;
0779   LogTrace("GsfElectronAlgo") << "E/p criteria are satisfied";
0780 
0781   // HoE cuts
0782   LogTrace("GsfElectronAlgo") << "HoE : " << ele.hcalOverEcal();
0783   double hoeCone = ele.hcalOverEcal();
0784   double hoeBc = ele.hcalOverEcalBc();
0785   const reco::CaloCluster& seedCluster = *(ele.superCluster()->seed());
0786   int detector = seedCluster.hitsAndFractions()[0].first.subdetId();
0787   bool HoEveto = false;
0788   double scle = ele.superCluster()->energy();
0789 
0790   if (detector == EcalBarrel)
0791     HoEveto = hoeCone * scle < cfg.maxHBarrelCone || hoeBc * scle < cfg.maxHBarrelBc ||
0792               hoeCone < cfg.maxHOverEBarrelCone || hoeBc < cfg.maxHOverEBarrelBc;
0793   else if (detector == EcalEndcap)
0794     HoEveto = hoeCone * scle < cfg.maxHEndcapsCone || hoeBc * scle < cfg.maxHEndcapsBc ||
0795               hoeCone < cfg.maxHOverEEndcapsCone || hoeBc < cfg.maxHOverEEndcapsBc;
0796 
0797   if (!HoEveto)
0798     return;
0799   LogTrace("GsfElectronAlgo") << "H/E criteria are satisfied";
0800 
0801   // delta eta criteria
0802   double deta = ele.deltaEtaSuperClusterTrackAtVtx();
0803   LogTrace("GsfElectronAlgo") << "delta eta : " << deta;
0804   if (ele.isEB() && (std::abs(deta) > cfg.maxDeltaEtaBarrel))
0805     return;
0806   if (ele.isEE() && (std::abs(deta) > cfg.maxDeltaEtaEndcaps))
0807     return;
0808   LogTrace("GsfElectronAlgo") << "Delta eta criteria are satisfied";
0809 
0810   // delta phi criteria
0811   double dphi = ele.deltaPhiSuperClusterTrackAtVtx();
0812   LogTrace("GsfElectronAlgo") << "delta phi : " << dphi;
0813   if (ele.isEB() && (std::abs(dphi) > cfg.maxDeltaPhiBarrel))
0814     return;
0815   if (ele.isEE() && (std::abs(dphi) > cfg.maxDeltaPhiEndcaps))
0816     return;
0817   LogTrace("GsfElectronAlgo") << "Delta phi criteria are satisfied";
0818 
0819   // sigma ieta ieta
0820   LogTrace("GsfElectronAlgo") << "sigma ieta ieta : " << ele.sigmaIetaIeta();
0821   if (ele.isEB() && (ele.sigmaIetaIeta() > cfg.maxSigmaIetaIetaBarrel))
0822     return;
0823   if (ele.isEE() && (ele.sigmaIetaIeta() > cfg.maxSigmaIetaIetaEndcaps))
0824     return;
0825   LogTrace("GsfElectronAlgo") << "Sigma ieta ieta criteria are satisfied";
0826 
0827   // fiducial
0828   if (!ele.isEB() && cfg.isBarrel)
0829     return;
0830   if (!ele.isEE() && cfg.isEndcaps)
0831     return;
0832   if (cfg.isFiducial &&
0833       (ele.isEBEEGap() || ele.isEBEtaGap() || ele.isEBPhiGap() || ele.isEERingGap() || ele.isEEDeeGap()))
0834     return;
0835   LogTrace("GsfElectronAlgo") << "Fiducial flags criteria are satisfied";
0836 
0837   // seed in TEC
0838   edm::RefToBase<TrajectorySeed> seed = ele.gsfTrack()->extra()->seedRef();
0839   ElectronSeedRef elseed = seed.castTo<ElectronSeedRef>();
0840   if (eg && !cfg_.cuts.seedFromTEC) {
0841     if (elseed.isNull()) {
0842       throw cms::Exception("GsfElectronAlgo|NotElectronSeed") << "The GsfTrack seed is not an ElectronSeed ?!";
0843     } else {
0844       if (elseed->subDet(1) == 6)
0845         return;
0846     }
0847   }
0848 
0849   // transverse impact parameter
0850   if (std::abs(ele.gsfTrack()->dxy(bs.position())) > cfg.maxTIP)
0851     return;
0852   LogTrace("GsfElectronAlgo") << "TIP criterion is satisfied";
0853 
0854   LogTrace("GsfElectronAlgo") << "All cut based criteria are satisfied";
0855   ele.setPassCutBasedPreselection(true);
0856 }
0857 
0858 void GsfElectronAlgo::createElectron(reco::GsfElectronCollection& electrons,
0859                                      ElectronData& electronData,
0860                                      EventData& eventData,
0861                                      CaloTopology const& topology,
0862                                      CaloGeometry const& geometry,
0863                                      MultiTrajectoryStateTransform const& mtsTransform,
0864                                      double magneticFieldInTesla,
0865                                      const GsfElectronAlgo::HeavyObjectCache* hoc,
0866                                      egamma::conv::TrackTableView ctfTable,
0867                                      egamma::conv::TrackTableView gsfTable,
0868                                      EcalPFRecHitThresholds const& thresholds,
0869                                      const HcalPFCuts* hcalCuts) {
0870   // charge ID
0871   int eleCharge;
0872   GsfElectron::ChargeInfo eleChargeInfo;
0873   electronData.computeCharge(eleCharge, eleChargeInfo);
0874 
0875   // electron basic cluster
0876   CaloClusterPtr elbcRef = electronData.getEleBasicCluster(mtsTransform);
0877 
0878   // Seed cluster
0879   const reco::CaloCluster& seedCluster = *(electronData.superClusterRef->seed());
0880 
0881   // seed Xtal
0882   // temporary, till CaloCluster->seed() is made available
0883   DetId seedXtalId = seedCluster.hitsAndFractions()[0].first;
0884 
0885   electronData.calculateMode();
0886 
0887   //====================================================
0888   // Candidate attributes
0889   //====================================================
0890 
0891   Candidate::LorentzVector momentum = electronData.calculateMomentum();
0892 
0893   //====================================================
0894   // Track-Cluster Matching
0895   //====================================================
0896 
0897   reco::GsfElectron::TrackClusterMatching tcMatching;
0898   tcMatching.electronCluster = elbcRef;
0899   tcMatching.eSuperClusterOverP =
0900       (electronData.vtxMom.mag() > 0) ? (electronData.superClusterRef->energy() / electronData.vtxMom.mag()) : (-1.);
0901   tcMatching.eSeedClusterOverP =
0902       (electronData.vtxMom.mag() > 0.) ? (seedCluster.energy() / electronData.vtxMom.mag()) : (-1);
0903   tcMatching.eSeedClusterOverPout =
0904       (electronData.seedMom.mag() > 0.) ? (seedCluster.energy() / electronData.seedMom.mag()) : (-1.);
0905   tcMatching.eEleClusterOverPout =
0906       (electronData.eleMom.mag() > 0.) ? (elbcRef->energy() / electronData.eleMom.mag()) : (-1.);
0907 
0908   EleRelPointPair scAtVtx(
0909       electronData.superClusterRef->position(), electronData.sclPos, eventData.beamspot->position());
0910   tcMatching.deltaEtaSuperClusterAtVtx = scAtVtx.dEta();
0911   tcMatching.deltaPhiSuperClusterAtVtx = scAtVtx.dPhi();
0912 
0913   EleRelPointPair seedAtCalo(seedCluster.position(), electronData.seedPos, eventData.beamspot->position());
0914   tcMatching.deltaEtaSeedClusterAtCalo = seedAtCalo.dEta();
0915   tcMatching.deltaPhiSeedClusterAtCalo = seedAtCalo.dPhi();
0916 
0917   EleRelPointPair ecAtCalo(elbcRef->position(), electronData.elePos, eventData.beamspot->position());
0918   tcMatching.deltaEtaEleClusterAtCalo = ecAtCalo.dEta();
0919   tcMatching.deltaPhiEleClusterAtCalo = ecAtCalo.dPhi();
0920 
0921   //=======================================================
0922   // Track extrapolations
0923   //=======================================================
0924 
0925   reco::GsfElectron::TrackExtrapolations tkExtra;
0926   ele_convert(electronData.vtxPos, tkExtra.positionAtVtx);
0927   ele_convert(electronData.sclPos, tkExtra.positionAtCalo);
0928   ele_convert(electronData.vtxMom, tkExtra.momentumAtVtx);
0929   ele_convert(electronData.sclMom, tkExtra.momentumAtCalo);
0930   ele_convert(electronData.seedMom, tkExtra.momentumOut);
0931   ele_convert(electronData.eleMom, tkExtra.momentumAtEleClus);
0932   ele_convert(electronData.vtxMomWithConstraint, tkExtra.momentumAtVtxWithConstraint);
0933 
0934   //=======================================================
0935   // Closest Ctf Track
0936   //=======================================================
0937 
0938   reco::GsfElectron::ClosestCtfTrack ctfInfo;
0939   ctfInfo.ctfTrack = electronData.ctfTrackRef;
0940   ctfInfo.shFracInnerHits = electronData.shFracInnerHits;
0941 
0942   //====================================================
0943   // FiducialFlags, using nextToBoundary definition of gaps
0944   //====================================================
0945 
0946   reco::GsfElectron::FiducialFlags fiducialFlags;
0947   int region = seedXtalId.det();
0948   int detector = seedXtalId.subdetId();
0949   double feta = std::abs(electronData.superClusterRef->position().eta());
0950   if (detector == EcalBarrel) {
0951     fiducialFlags.isEB = true;
0952     EBDetId ebdetid(seedXtalId);
0953     if (EBDetId::isNextToEtaBoundary(ebdetid)) {
0954       if (ebdetid.ietaAbs() == 85) {
0955         fiducialFlags.isEBEEGap = true;
0956       } else {
0957         fiducialFlags.isEBEtaGap = true;
0958       }
0959     }
0960     if (EBDetId::isNextToPhiBoundary(ebdetid)) {
0961       fiducialFlags.isEBPhiGap = true;
0962     }
0963   } else if (detector == EcalEndcap) {
0964     fiducialFlags.isEE = true;
0965     EEDetId eedetid(seedXtalId);
0966     if (EEDetId::isNextToRingBoundary(eedetid)) {
0967       if (std::abs(feta) < 2.) {
0968         fiducialFlags.isEBEEGap = true;
0969       } else {
0970         fiducialFlags.isEERingGap = true;
0971       }
0972     }
0973     if (EEDetId::isNextToDBoundary(eedetid)) {
0974       fiducialFlags.isEEDeeGap = true;
0975     }
0976   } else if (EcalTools::isHGCalDet((DetId::Detector)region)) {
0977     fiducialFlags.isEE = true;
0978     //HGCalDetId eeDetid(seedXtalId);
0979     // fill in fiducial information when we know how to use it...
0980   } else {
0981     throw cms::Exception("GsfElectronAlgo|UnknownXtalRegion")
0982         << "createElectron(): do not know if it is a barrel or endcap seed cluster !!!!";
0983   }
0984 
0985   //====================================================
0986   // SaturationInfo
0987   //====================================================
0988 
0989   auto saturationInfo = calculateSaturationInfo(electronData.superClusterRef, eventData);
0990 
0991   //====================================================
0992   // ShowerShape
0993   //====================================================
0994 
0995   reco::GsfElectron::ShowerShape showerShape;
0996   reco::GsfElectron::ShowerShape full5x5_showerShape;
0997   if (!EcalTools::isHGCalDet((DetId::Detector)region)) {
0998     showerShape = calculateShowerShape<false>(electronData.superClusterRef,
0999                                               hcalHelperCone_,
1000                                               hcalHelperBc_,
1001                                               eventData,
1002                                               topology,
1003                                               geometry,
1004                                               thresholds,
1005                                               hcalCuts);
1006     full5x5_showerShape = calculateShowerShape<true>(electronData.superClusterRef,
1007                                                      hcalHelperCone_,
1008                                                      hcalHelperBc_,
1009                                                      eventData,
1010                                                      topology,
1011                                                      geometry,
1012                                                      thresholds,
1013                                                      hcalCuts);
1014   }
1015 
1016   //====================================================
1017   // ConversionRejection
1018   //====================================================
1019 
1020   edm::Handle<reco::TrackCollection> ctfTracks = eventData.originalCtfTracks;
1021   if (!ctfTracks.isValid()) {
1022     ctfTracks = eventData.currentCtfTracks;
1023   }
1024 
1025   {
1026     //get the references to the gsf and ctf tracks that are made
1027     //by the electron
1028     const reco::TrackRef el_ctftrack = electronData.coreRef->ctfTrack();
1029     const reco::GsfTrackRef& el_gsftrack = electronData.coreRef->gsfTrack();
1030 
1031     //protect against the wrong collection being passed to the function
1032     if (el_ctftrack.isNonnull() && el_ctftrack.id() != ctfTracks.id())
1033       throw cms::Exception("ConversionFinderError")
1034           << "ProductID of ctf track collection does not match ProductID of electron's CTF track! \n";
1035     if (el_gsftrack.isNonnull() && el_gsftrack.id() != eventData.originalGsfTracks.id())
1036       throw cms::Exception("ConversionFinderError")
1037           << "ProductID of gsf track collection does not match ProductID of electron's GSF track! \n";
1038   }
1039 
1040   // values of conversionInfo.flag
1041   // -9999 : Partner track was not found
1042   // 0     : Partner track found in the CTF collection using
1043   // 1     : Partner track found in the CTF collection using
1044   // 2     : Partner track found in the GSF collection using
1045   // 3     : Partner track found in the GSF collection using the electron's GSF track
1046   auto conversionInfo = egamma::conv::findConversion(*electronData.coreRef, ctfTable, gsfTable, magneticFieldInTesla);
1047 
1048   reco::GsfElectron::ConversionRejection conversionVars;
1049   conversionVars.flags = conversionInfo.flag;
1050   conversionVars.dist = conversionInfo.dist;
1051   conversionVars.dcot = conversionInfo.dcot;
1052   conversionVars.radius = conversionInfo.radiusOfConversion;
1053   if (cfg_.strategy.fillConvVtxFitProb) {
1054     //this is an intentionally bugged version which ignores the GsfTrack
1055     //this is a bug which was introduced in reduced e/gamma where the GsfTrack gets
1056     //relinked to a new collection which means it can no longer match the conversion
1057     //as it matches based on product/id
1058     //we keep this defination for the MVAs
1059     const auto matchedConv = ConversionTools::matchedConversion(
1060         electronData.coreRef->ctfTrack(), *eventData.conversions, eventData.beamspot->position(), 2.0, 1e-6, 0);
1061     conversionVars.vtxFitProb = ConversionTools::getVtxFitProb(matchedConv);
1062   }
1063   if (conversionInfo.conversionPartnerCtfTkIdx) {
1064     conversionVars.partner = TrackBaseRef(reco::TrackRef(ctfTracks, conversionInfo.conversionPartnerCtfTkIdx.value()));
1065   } else if (conversionInfo.conversionPartnerGsfTkIdx) {
1066     conversionVars.partner =
1067         TrackBaseRef(reco::GsfTrackRef(eventData.originalGsfTracks, conversionInfo.conversionPartnerGsfTkIdx.value()));
1068   }
1069 
1070   //====================================================
1071   // Go !
1072   //====================================================
1073 
1074   electrons.emplace_back(eleCharge,
1075                          eleChargeInfo,
1076                          electronData.coreRef,
1077                          tcMatching,
1078                          tkExtra,
1079                          ctfInfo,
1080                          fiducialFlags,
1081                          showerShape,
1082                          full5x5_showerShape,
1083                          conversionVars,
1084                          saturationInfo);
1085   auto& ele = electrons.back();
1086   // Will be overwritten later in the case of the regression
1087   ele.setCorrectedEcalEnergyError(egamma::ecalClusterEnergyUncertaintyElectronSpecific(*(ele.superCluster())));
1088   ele.setP4(GsfElectron::P4_FROM_SUPER_CLUSTER, momentum, 0, true);
1089   ele.setMass(0.0);
1090 
1091   //====================================================
1092   // brems fractions
1093   //====================================================
1094 
1095   if (electronData.innMom.mag() > 0.) {
1096     ele.setTrackFbrem((electronData.innMom.mag() - electronData.outMom.mag()) / electronData.innMom.mag());
1097   }
1098 
1099   // the supercluster is the refined one The seed is not necessarily the first cluster
1100   // hence the use of the electronCluster
1101   SuperClusterRef sc = ele.superCluster();
1102   if (!(sc.isNull())) {
1103     CaloClusterPtr cl = ele.electronCluster();
1104     if (sc->clustersSize() > 1) {
1105       float pf_fbrem = (sc->energy() - cl->energy()) / sc->energy();
1106       ele.setSuperClusterFbrem(pf_fbrem);
1107     } else {
1108       ele.setSuperClusterFbrem(0);
1109     }
1110   }
1111 
1112   //====================================================
1113   // classification and corrections
1114   //====================================================
1115   // classification
1116   const auto elClass = egamma::classifyElectron(ele);
1117   ele.setClassification(elClass);
1118 
1119   bool unexpectedClassification = elClass == GsfElectron::UNKNOWN || elClass > GsfElectron::GAP;
1120   if (unexpectedClassification) {
1121     edm::LogWarning("GsfElectronAlgo") << "unexpected classification";
1122   }
1123 
1124   // ecal energy
1125   if (cfg_.strategy.useEcalRegression)  // new
1126   {
1127     regHelper_.applyEcalRegression(ele, *eventData.vertices, *eventData.barrelRecHits, *eventData.endcapRecHits);
1128   } else  // original implementation
1129   {
1130     if (!EcalTools::isHGCalDet((DetId::Detector)region)) {
1131       if (ele.core()->ecalDrivenSeed()) {
1132         if (cfg_.strategy.ecalDrivenEcalEnergyFromClassBasedParameterization && !unexpectedClassification) {
1133           if (ele.isEcalEnergyCorrected()) {
1134             edm::LogWarning("ElectronEnergyCorrector::classBasedElectronEnergy") << "already done";
1135           } else {
1136             ele.setCorrectedEcalEnergy(
1137                 egamma::classBasedElectronEnergy(ele, *eventData.beamspot, *crackCorrectionFunction_));
1138           }
1139         }
1140         if (cfg_.strategy.ecalDrivenEcalErrorFromClassBasedParameterization) {
1141           ele.setCorrectedEcalEnergyError(egamma::classBasedElectronEnergyUncertainty(ele));
1142         }
1143       } else {
1144         if (cfg_.strategy.pureTrackerDrivenEcalErrorFromSimpleParameterization) {
1145           ele.setCorrectedEcalEnergyError(egamma::simpleElectronEnergyUncertainty(ele));
1146         }
1147       }
1148     }
1149   }
1150 
1151   // momentum
1152   // Keep the default correction running first. The track momentum error is computed in there
1153   if (cfg_.strategy.useDefaultEnergyCorrection && ele.core()->ecalDrivenSeed() && !unexpectedClassification) {
1154     if (ele.p4Error(reco::GsfElectron::P4_COMBINATION) != 999.) {
1155       edm::LogWarning("ElectronMomentumCorrector::correct") << "already done";
1156     } else {
1157       auto p = egamma::correctElectronMomentum(ele, electronData.vtxTSOS);
1158       ele.correctMomentum(p.momentum, p.trackError, p.finalError);
1159     }
1160   }
1161   if (cfg_.strategy.useCombinationRegression)  // new
1162   {
1163     regHelper_.applyCombinationRegression(ele);
1164   }
1165 
1166   //====================================================
1167   // now isolation variables
1168   //====================================================
1169   reco::GsfElectron::IsolationVariables dr03, dr04;
1170   dr03.tkSumPt = eventData.tkIsol03Calc(*ele.gsfTrack()).ptSum;
1171   dr04.tkSumPt = eventData.tkIsol04Calc(*ele.gsfTrack()).ptSum;
1172   dr03.tkSumPtHEEP = eventData.tkIsolHEEP03Calc(*ele.gsfTrack()).ptSum;
1173   dr04.tkSumPtHEEP = eventData.tkIsolHEEP04Calc(*ele.gsfTrack()).ptSum;
1174 
1175   if (!EcalTools::isHGCalDet((DetId::Detector)region)) {
1176     for (uint id = 0; id < dr03.hcalRecHitSumEt.size(); ++id) {
1177       dr03.hcalRecHitSumEt[id] = eventData.hadIsolation03.getHcalEtSum(&ele, id + 1, hcalCuts);
1178       dr03.hcalRecHitSumEtBc[id] = eventData.hadIsolation03Bc.getHcalEtSumBc(&ele, id + 1, hcalCuts);
1179 
1180       dr04.hcalRecHitSumEt[id] = eventData.hadIsolation04.getHcalEtSum(&ele, id + 1, hcalCuts);
1181       dr04.hcalRecHitSumEtBc[id] = eventData.hadIsolation04Bc.getHcalEtSumBc(&ele, id + 1, hcalCuts);
1182     }
1183 
1184     dr03.ecalRecHitSumEt = eventData.ecalBarrelIsol03.getEtSum(&ele, thresholds);
1185     dr03.ecalRecHitSumEt += eventData.ecalEndcapIsol03.getEtSum(&ele, thresholds);
1186 
1187     dr04.ecalRecHitSumEt = eventData.ecalBarrelIsol04.getEtSum(&ele, thresholds);
1188     dr04.ecalRecHitSumEt += eventData.ecalEndcapIsol04.getEtSum(&ele, thresholds);
1189   }
1190 
1191   dr03.pre7DepthHcal = false;
1192   dr04.pre7DepthHcal = false;
1193 
1194   ele.setIsolation03(dr03);
1195   ele.setIsolation04(dr04);
1196 
1197   //====================================================
1198   // PFclusters based ISO !
1199   // Added in CMSSW_12_0_1 at this stage to be able to use them in PF ID DNN
1200   //====================================================
1201   if (cfg_.strategy.computePfClusterIso) {
1202     reco::GsfElectron::PflowIsolationVariables isoVariables;
1203     isoVariables.sumEcalClusterEt = ecalisoAlgo_->getSum(ele, eventData.ecalClustersHandle);
1204     isoVariables.sumHcalClusterEt = hcalisoAlgo_->getSum(ele, eventData.hcalClustersHandle);
1205     // Other Pfiso variables are initialized at 0 and not used
1206     ele.setPfIsolationVariables(isoVariables);
1207   }
1208 
1209   //====================================================
1210   // preselection flag
1211   //====================================================
1212 
1213   setCutBasedPreselectionFlag(ele, *eventData.beamspot);
1214   //setting mva flag, currently GedGsfElectron and GsfElectron pre-selection flags have desynced
1215   //this is for GedGsfElectrons, GsfElectrons (ie old pre 7X std reco) resets this later on
1216   //in the function "addPfInfo"
1217   //yes this is awful, we'll fix it once we work out how to...
1218   float mvaValue = hoc->sElectronMVAEstimator->mva(ele, *(eventData.vertices));
1219   ele.setPassMvaPreselection(mvaValue > cfg_.strategy.PreSelectMVA);
1220 
1221   //====================================================
1222   // Pixel match variables
1223   //====================================================
1224   setPixelMatchInfomation(ele);
1225 
1226   LogTrace("GsfElectronAlgo") << "Constructed new electron with energy  " << ele.p4().e();
1227 }
1228 
1229 // Pixel match variables
1230 void GsfElectronAlgo::setPixelMatchInfomation(reco::GsfElectron& ele) const {
1231   int sd1 = 0;
1232   int sd2 = 0;
1233   float dPhi1 = 0;
1234   float dPhi2 = 0;
1235   float dRz1 = 0;
1236   float dRz2 = 0;
1237   edm::RefToBase<TrajectorySeed> seed = ele.gsfTrack()->extra()->seedRef();
1238   ElectronSeedRef elseed = seed.castTo<ElectronSeedRef>();
1239   if (seed.isNull()) {
1240   } else {
1241     if (elseed.isNull()) {
1242     } else {
1243       sd1 = elseed->subDet(0);
1244       sd2 = elseed->subDet(1);
1245       dPhi1 = (ele.charge() > 0) ? elseed->dPhiPos(0) : elseed->dPhiNeg(0);
1246       dPhi2 = (ele.charge() > 0) ? elseed->dPhiPos(1) : elseed->dPhiNeg(1);
1247       dRz1 = (ele.charge() > 0) ? elseed->dRZPos(0) : elseed->dRZNeg(0);
1248       dRz2 = (ele.charge() > 0) ? elseed->dRZPos(1) : elseed->dRZNeg(1);
1249     }
1250   }
1251   ele.setPixelMatchSubdetectors(sd1, sd2);
1252   ele.setPixelMatchDPhi1(dPhi1);
1253   ele.setPixelMatchDPhi2(dPhi2);
1254   ele.setPixelMatchDRz1(dRz1);
1255   ele.setPixelMatchDRz2(dRz2);
1256 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!