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
0041 SoftElectronMVAEstimator::Configuration sconfig;
0042 sconfig.vweightsfiles = conf.getParameter<std::vector<std::string>>("SoftElecMVAFilesString");
0043 sElectronMVAEstimator = std::make_unique<SoftElectronMVAEstimator>(sconfig);
0044
0045 ElectronMVAEstimator::Configuration iconfig;
0046 iconfig.vweightsfiles = conf.getParameter<std::vector<std::string>>("ElecMVAFilesString");
0047 iElectronMVAEstimator = std::make_unique<ElectronMVAEstimator>(iconfig);
0048
0049
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
0066
0067
0068 struct GsfElectronAlgo::EventData {
0069
0070 void retreiveOriginalTrackCollections(const reco::TrackRef&, const reco::GsfTrackRef&);
0071
0072
0073 edm::Event const* event;
0074 const reco::BeamSpot* beamspot;
0075
0076
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
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
0103 edm::Handle<reco::PFClusterCollection> ecalClustersHandle;
0104 std::vector<edm::Handle<reco::PFClusterCollection>> hcalClustersHandle;
0105 };
0106
0107
0108
0109
0110
0111 struct GsfElectronAlgo::ElectronData {
0112
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
0121 ElectronData(const reco::GsfElectronCoreRef& core, const reco::BeamSpot& bs);
0122
0123
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
0131 TrajectoryStateOnSurface innTSOS;
0132 TrajectoryStateOnSurface outTSOS;
0133 TrajectoryStateOnSurface vtxTSOS;
0134 TrajectoryStateOnSurface sclTSOS;
0135 TrajectoryStateOnSurface seedTSOS;
0136 TrajectoryStateOnSurface eleTSOS;
0137 TrajectoryStateOnSurface constrainedVtxTSOS;
0138
0139
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
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
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
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
0217 innTSOS = mtsTransform.innerStateOnSurface(*gsfTrackRef);
0218 if (!innTSOS.isValid())
0219 return false;
0220
0221
0222
0223 GlobalPoint bsPos;
0224 ele_convert(beamSpot.position(), bsPos);
0225 vtxTSOS = mtsTransform.extrapolatedState(innTSOS, bsPos);
0226 if (!vtxTSOS.isValid())
0227 vtxTSOS = innTSOS;
0228
0229
0230 outTSOS = mtsTransform.outerStateOnSurface(*gsfTrackRef);
0231 if (!outTSOS.isValid())
0232 return false;
0233
0234
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
0243 sclTSOS = mtsTransform.extrapolatedState(
0244 innTSOS, GlobalPoint(superClusterRef->x(), superClusterRef->y(), superClusterRef->z()));
0245 if (!sclTSOS.isValid())
0246 sclTSOS = outTSOS;
0247
0248
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
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
0339
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
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
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
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
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
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
0692 const GsfElectronCoreRef coreRef = edm::Ref<GsfElectronCoreCollection>(eventData.coreElectrons, i);
0693
0694
0695 if (coreRef->superCluster().isNull())
0696 continue;
0697
0698
0699 ElectronData electronData(coreRef, *eventData.beamspot);
0700
0701
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 }
0731 return electrons;
0732 }
0733
0734 void GsfElectronAlgo::setCutBasedPreselectionFlag(GsfElectron& ele, const reco::BeamSpot& bs) const {
0735
0736 ele.setPassCutBasedPreselection(false);
0737
0738
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
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
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
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
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
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
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
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
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
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
0864 int eleCharge;
0865 GsfElectron::ChargeInfo eleChargeInfo;
0866 electronData.computeCharge(eleCharge, eleChargeInfo);
0867
0868
0869 CaloClusterPtr elbcRef = electronData.getEleBasicCluster(mtsTransform);
0870
0871
0872 const reco::CaloCluster& seedCluster = *(electronData.superClusterRef->seed());
0873
0874
0875
0876 DetId seedXtalId = seedCluster.hitsAndFractions()[0].first;
0877
0878 electronData.calculateMode();
0879
0880
0881
0882
0883
0884 Candidate::LorentzVector momentum = electronData.calculateMomentum();
0885
0886
0887
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
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
0929
0930
0931 reco::GsfElectron::ClosestCtfTrack ctfInfo;
0932 ctfInfo.ctfTrack = electronData.ctfTrackRef;
0933 ctfInfo.shFracInnerHits = electronData.shFracInnerHits;
0934
0935
0936
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
0972
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
0980
0981
0982 auto saturationInfo = calculateSaturationInfo(electronData.superClusterRef, eventData);
0983
0984
0985
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
0999
1000
1001 edm::Handle<reco::TrackCollection> ctfTracks = eventData.originalCtfTracks;
1002 if (!ctfTracks.isValid()) {
1003 ctfTracks = eventData.currentCtfTracks;
1004 }
1005
1006 {
1007
1008
1009 const reco::TrackRef el_ctftrack = electronData.coreRef->ctfTrack();
1010 const reco::GsfTrackRef& el_gsftrack = electronData.coreRef->gsfTrack();
1011
1012
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
1022
1023
1024
1025
1026
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
1036
1037
1038
1039
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
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
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
1074
1075
1076 if (electronData.innMom.mag() > 0.) {
1077 ele.setTrackFbrem((electronData.innMom.mag() - electronData.outMom.mag()) / electronData.innMom.mag());
1078 }
1079
1080
1081
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
1095
1096
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
1106 if (cfg_.strategy.useEcalRegression)
1107 {
1108 regHelper_.applyEcalRegression(ele, *eventData.vertices, *eventData.barrelRecHits, *eventData.endcapRecHits);
1109 } else
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
1133
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)
1143 {
1144 regHelper_.applyCombinationRegression(ele);
1145 }
1146
1147
1148
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
1180
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
1187 ele.setPfIsolationVariables(isoVariables);
1188 }
1189
1190
1191
1192
1193
1194 setCutBasedPreselectionFlag(ele, *eventData.beamspot);
1195
1196
1197
1198
1199 float mvaValue = hoc->sElectronMVAEstimator->mva(ele, *(eventData.vertices));
1200 ele.setPassMvaPreselection(mvaValue > cfg_.strategy.PreSelectMVA);
1201
1202
1203
1204
1205 setPixelMatchInfomation(ele);
1206
1207 LogTrace("GsfElectronAlgo") << "Constructed new electron with energy " << ele.p4().e();
1208 }
1209
1210
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 }