File indexing completed on 2024-10-08 05:11:57
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
0044 SoftElectronMVAEstimator::Configuration sconfig;
0045 sconfig.vweightsfiles = conf.getParameter<std::vector<std::string>>("SoftElecMVAFilesString");
0046 sElectronMVAEstimator = std::make_unique<SoftElectronMVAEstimator>(sconfig);
0047
0048 ElectronMVAEstimator::Configuration iconfig;
0049 iconfig.vweightsfiles = conf.getParameter<std::vector<std::string>>("ElecMVAFilesString");
0050 iElectronMVAEstimator = std::make_unique<ElectronMVAEstimator>(iconfig);
0051
0052
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
0069
0070
0071 struct GsfElectronAlgo::EventData {
0072
0073 void retreiveOriginalTrackCollections(const reco::TrackRef&, const reco::GsfTrackRef&);
0074
0075
0076 edm::Event const* event;
0077 const reco::BeamSpot* beamspot;
0078
0079
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
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
0106 edm::Handle<reco::PFClusterCollection> ecalClustersHandle;
0107 std::vector<edm::Handle<reco::PFClusterCollection>> hcalClustersHandle;
0108 };
0109
0110
0111
0112
0113
0114 struct GsfElectronAlgo::ElectronData {
0115
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
0124 ElectronData(const reco::GsfElectronCoreRef& core, const reco::BeamSpot& bs);
0125
0126
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
0134 TrajectoryStateOnSurface innTSOS;
0135 TrajectoryStateOnSurface outTSOS;
0136 TrajectoryStateOnSurface vtxTSOS;
0137 TrajectoryStateOnSurface sclTSOS;
0138 TrajectoryStateOnSurface seedTSOS;
0139 TrajectoryStateOnSurface eleTSOS;
0140 TrajectoryStateOnSurface constrainedVtxTSOS;
0141
0142
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
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
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
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
0220 innTSOS = mtsTransform.innerStateOnSurface(*gsfTrackRef);
0221 if (!innTSOS.isValid())
0222 return false;
0223
0224
0225
0226 GlobalPoint bsPos;
0227 ele_convert(beamSpot.position(), bsPos);
0228 vtxTSOS = mtsTransform.extrapolatedState(innTSOS, bsPos);
0229 if (!vtxTSOS.isValid())
0230 vtxTSOS = innTSOS;
0231
0232
0233 outTSOS = mtsTransform.outerStateOnSurface(*gsfTrackRef);
0234 if (!outTSOS.isValid())
0235 return false;
0236
0237
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
0246 sclTSOS = mtsTransform.extrapolatedState(
0247 innTSOS, GlobalPoint(superClusterRef->x(), superClusterRef->y(), superClusterRef->z()));
0248 if (!sclTSOS.isValid())
0249 sclTSOS = outTSOS;
0250
0251
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
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
0343
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
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
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
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
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
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
0697 const GsfElectronCoreRef coreRef = edm::Ref<GsfElectronCoreCollection>(eventData.coreElectrons, i);
0698
0699
0700 if (coreRef->superCluster().isNull())
0701 continue;
0702
0703
0704 ElectronData electronData(coreRef, *eventData.beamspot);
0705
0706
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 }
0737 return electrons;
0738 }
0739
0740 void GsfElectronAlgo::setCutBasedPreselectionFlag(GsfElectron& ele, const reco::BeamSpot& bs) const {
0741
0742 ele.setPassCutBasedPreselection(false);
0743
0744
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
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
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
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
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
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
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
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
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
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
0871 int eleCharge;
0872 GsfElectron::ChargeInfo eleChargeInfo;
0873 electronData.computeCharge(eleCharge, eleChargeInfo);
0874
0875
0876 CaloClusterPtr elbcRef = electronData.getEleBasicCluster(mtsTransform);
0877
0878
0879 const reco::CaloCluster& seedCluster = *(electronData.superClusterRef->seed());
0880
0881
0882
0883 DetId seedXtalId = seedCluster.hitsAndFractions()[0].first;
0884
0885 electronData.calculateMode();
0886
0887
0888
0889
0890
0891 Candidate::LorentzVector momentum = electronData.calculateMomentum();
0892
0893
0894
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
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
0936
0937
0938 reco::GsfElectron::ClosestCtfTrack ctfInfo;
0939 ctfInfo.ctfTrack = electronData.ctfTrackRef;
0940 ctfInfo.shFracInnerHits = electronData.shFracInnerHits;
0941
0942
0943
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
0979
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
0987
0988
0989 auto saturationInfo = calculateSaturationInfo(electronData.superClusterRef, eventData);
0990
0991
0992
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
1018
1019
1020 edm::Handle<reco::TrackCollection> ctfTracks = eventData.originalCtfTracks;
1021 if (!ctfTracks.isValid()) {
1022 ctfTracks = eventData.currentCtfTracks;
1023 }
1024
1025 {
1026
1027
1028 const reco::TrackRef el_ctftrack = electronData.coreRef->ctfTrack();
1029 const reco::GsfTrackRef& el_gsftrack = electronData.coreRef->gsfTrack();
1030
1031
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
1041
1042
1043
1044
1045
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
1055
1056
1057
1058
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
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
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
1093
1094
1095 if (electronData.innMom.mag() > 0.) {
1096 ele.setTrackFbrem((electronData.innMom.mag() - electronData.outMom.mag()) / electronData.innMom.mag());
1097 }
1098
1099
1100
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
1114
1115
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
1125 if (cfg_.strategy.useEcalRegression)
1126 {
1127 regHelper_.applyEcalRegression(ele, *eventData.vertices, *eventData.barrelRecHits, *eventData.endcapRecHits);
1128 } else
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
1152
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)
1162 {
1163 regHelper_.applyCombinationRegression(ele);
1164 }
1165
1166
1167
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
1199
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
1206 ele.setPfIsolationVariables(isoVariables);
1207 }
1208
1209
1210
1211
1212
1213 setCutBasedPreselectionFlag(ele, *eventData.beamspot);
1214
1215
1216
1217
1218 float mvaValue = hoc->sElectronMVAEstimator->mva(ele, *(eventData.vertices));
1219 ele.setPassMvaPreselection(mvaValue > cfg_.strategy.PreSelectMVA);
1220
1221
1222
1223
1224 setPixelMatchInfomation(ele);
1225
1226 LogTrace("GsfElectronAlgo") << "Constructed new electron with energy " << ele.p4().e();
1227 }
1228
1229
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 }