Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-05-15 22:20:47 

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