File indexing completed on 2023-10-25 09:59:09
0001 #include "RecoEgamma/EgammaElectronProducers/interface/LowPtGsfElectronFeatures.h"
0002 #include "CommonTools/BaseParticlePropagator/interface/BaseParticlePropagator.h"
0003 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
0004 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
0005 #include "DataFormats/ParticleFlowReco/interface/PFClusterFwd.h"
0006 #include "DataFormats/PatCandidates/interface/Electron.h"
0007 #include "DataFormats/TrackReco/interface/Track.h"
0008 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0009 #include "TVector3.h"
0010 #include <algorithm>
0011 #include <iostream>
0012
0013 namespace lowptgsfeleseed {
0014
0015
0016
0017 std::vector<float> features(const reco::PreId& ecal,
0018 const reco::PreId& hcal,
0019 double rho,
0020 const reco::BeamSpot& spot,
0021 noZS::EcalClusterLazyTools& tools) {
0022 float trk_pt_ = -1.;
0023 float trk_eta_ = -1.;
0024 float trk_phi_ = -1.;
0025 float trk_p_ = -1.;
0026 float trk_nhits_ = -1.;
0027 float trk_high_quality_ = -1.;
0028 float trk_chi2red_ = -1.;
0029 float rho_ = -1.;
0030 float ktf_ecal_cluster_e_ = -1.;
0031 float ktf_ecal_cluster_deta_ = -42.;
0032 float ktf_ecal_cluster_dphi_ = -42.;
0033 float ktf_ecal_cluster_e3x3_ = -1.;
0034 float ktf_ecal_cluster_e5x5_ = -1.;
0035 float ktf_ecal_cluster_covEtaEta_ = -42.;
0036 float ktf_ecal_cluster_covEtaPhi_ = -42.;
0037 float ktf_ecal_cluster_covPhiPhi_ = -42.;
0038 float ktf_ecal_cluster_r9_ = -0.1;
0039 float ktf_ecal_cluster_circularity_ = -0.1;
0040 float ktf_hcal_cluster_e_ = -1.;
0041 float ktf_hcal_cluster_deta_ = -42.;
0042 float ktf_hcal_cluster_dphi_ = -42.;
0043 float preid_gsf_dpt_ = -1.;
0044 float preid_trk_gsf_chiratio_ = -1.;
0045 float preid_gsf_chi2red_ = -1.;
0046 float trk_dxy_sig_ = -1.;
0047
0048
0049 const auto& trk = ecal.trackRef();
0050 if (trk.isNonnull()) {
0051 trk_pt_ = trk->pt();
0052 trk_eta_ = trk->eta();
0053 trk_phi_ = trk->phi();
0054 trk_p_ = trk->p();
0055 trk_nhits_ = static_cast<float>(trk->found());
0056 trk_high_quality_ = static_cast<float>(trk->quality(reco::TrackBase::qualityByName("highPurity")));
0057 trk_chi2red_ = trk->normalizedChi2();
0058 if (trk->dxy(spot) > 0.) {
0059 trk_dxy_sig_ = trk->dxyError() / trk->dxy(spot);
0060 }
0061 ktf_ecal_cluster_dphi_ *= trk->charge();
0062 }
0063
0064
0065 rho_ = static_cast<float>(rho);
0066
0067
0068 const auto& ecal_clu = ecal.clusterRef();
0069 if (ecal_clu.isNonnull()) {
0070 ktf_ecal_cluster_e_ = ecal_clu->energy();
0071 ktf_ecal_cluster_deta_ = ecal.geomMatching()[0];
0072 ktf_ecal_cluster_dphi_ = ecal.geomMatching()[1];
0073 ktf_ecal_cluster_e3x3_ = tools.e3x3(*ecal_clu);
0074 ktf_ecal_cluster_e5x5_ = tools.e5x5(*ecal_clu);
0075 const auto& covs = tools.localCovariances(*ecal_clu);
0076 ktf_ecal_cluster_covEtaEta_ = covs[0];
0077 ktf_ecal_cluster_covEtaPhi_ = covs[1];
0078 ktf_ecal_cluster_covPhiPhi_ = covs[2];
0079 if (ktf_ecal_cluster_e_ > 0.) {
0080 ktf_ecal_cluster_r9_ = ktf_ecal_cluster_e3x3_ / ktf_ecal_cluster_e_;
0081 }
0082 if (ktf_ecal_cluster_e5x5_ > 0.) {
0083 ktf_ecal_cluster_circularity_ = 1. - tools.e1x5(*ecal_clu) / ktf_ecal_cluster_e5x5_;
0084 } else {
0085 ktf_ecal_cluster_circularity_ = -0.1;
0086 }
0087 }
0088
0089
0090 const auto& hcal_clu = hcal.clusterRef();
0091 if (hcal_clu.isNonnull()) {
0092 ktf_hcal_cluster_e_ = hcal_clu->energy();
0093 ktf_hcal_cluster_deta_ = hcal.geomMatching()[0];
0094 ktf_hcal_cluster_dphi_ = hcal.geomMatching()[1];
0095 }
0096
0097
0098 preid_gsf_dpt_ = ecal.dpt();
0099 preid_trk_gsf_chiratio_ = ecal.chi2Ratio();
0100 preid_gsf_chi2red_ = ecal.gsfChi2();
0101
0102
0103 std::vector<float> output = {trk_pt_,
0104 trk_eta_,
0105 trk_phi_,
0106 trk_p_,
0107 trk_nhits_,
0108 trk_high_quality_,
0109 trk_chi2red_,
0110 rho_,
0111 ktf_ecal_cluster_e_,
0112 ktf_ecal_cluster_deta_,
0113 ktf_ecal_cluster_dphi_,
0114 ktf_ecal_cluster_e3x3_,
0115 ktf_ecal_cluster_e5x5_,
0116 ktf_ecal_cluster_covEtaEta_,
0117 ktf_ecal_cluster_covEtaPhi_,
0118 ktf_ecal_cluster_covPhiPhi_,
0119 ktf_ecal_cluster_r9_,
0120 ktf_ecal_cluster_circularity_,
0121 ktf_hcal_cluster_e_,
0122 ktf_hcal_cluster_deta_,
0123 ktf_hcal_cluster_dphi_,
0124 preid_gsf_dpt_,
0125 preid_trk_gsf_chiratio_,
0126 preid_gsf_chi2red_,
0127 trk_dxy_sig_};
0128 return output;
0129 };
0130
0131 }
0132
0133 namespace lowptgsfeleid {
0134
0135 std::vector<float> features_V1(
0136 reco::GsfElectron const& ele, float rho, float unbiased, float field_z, const reco::Track* trk) {
0137 float eid_rho = -999.;
0138 float eid_sc_eta = -999.;
0139 float eid_shape_full5x5_r9 = -999.;
0140 float eid_sc_etaWidth = -999.;
0141 float eid_sc_phiWidth = -999.;
0142 float eid_shape_full5x5_HoverE = -999.;
0143 float eid_trk_nhits = -999.;
0144 float eid_trk_chi2red = -999.;
0145 float eid_gsf_chi2red = -999.;
0146 float eid_brem_frac = -999.;
0147 float eid_gsf_nhits = -999.;
0148 float eid_match_SC_EoverP = -999.;
0149 float eid_match_eclu_EoverP = -999.;
0150 float eid_match_SC_dEta = -999.;
0151 float eid_match_SC_dPhi = -999.;
0152 float eid_match_seed_dEta = -999.;
0153 float eid_sc_E = -999.;
0154 float eid_trk_p = -999.;
0155 float gsf_mode_p = -999.;
0156 float core_shFracHits = -999.;
0157 float gsf_bdtout1 = -999.;
0158 float gsf_dr = -999.;
0159 float trk_dr = -999.;
0160 float sc_Nclus = -999.;
0161 float sc_clus1_nxtal = -999.;
0162 float sc_clus1_dphi = -999.;
0163 float sc_clus2_dphi = -999.;
0164 float sc_clus1_deta = -999.;
0165 float sc_clus2_deta = -999.;
0166 float sc_clus1_E = -999.;
0167 float sc_clus2_E = -999.;
0168 float sc_clus1_E_ov_p = -999.;
0169 float sc_clus2_E_ov_p = -999.;
0170
0171
0172 if (trk != nullptr || (ele.core().isNonnull() && ele.closestCtfTrackRef().isNonnull())) {
0173 const reco::Track* tk = trk ? trk : &*ele.closestCtfTrackRef();
0174 eid_trk_p = tk->p();
0175 eid_trk_nhits = tk->found();
0176 eid_trk_chi2red = tk->normalizedChi2();
0177 trk_dr = reco::deltaR(*tk, ele);
0178 }
0179
0180
0181 if (ele.core().isNonnull()) {
0182 reco::GsfTrackRef gsf = ele.gsfTrack();
0183 if (gsf.isNonnull()) {
0184 gsf_mode_p = gsf->pMode();
0185 eid_gsf_nhits = (float)gsf->found();
0186 eid_gsf_chi2red = gsf->normalizedChi2();
0187 TVector3 gsfTV3(0, 0, 0);
0188 gsfTV3.SetPtEtaPhi(gsf->ptMode(), gsf->etaMode(), gsf->phiMode());
0189 TVector3 eleTV3(0, 0, 0);
0190 eleTV3.SetPtEtaPhi(ele.pt(), ele.eta(), ele.phi());
0191 gsf_dr = eleTV3.DeltaR(gsfTV3);
0192 }
0193 }
0194
0195
0196 if (ele.core().isNonnull()) {
0197 reco::SuperClusterRef sc = ele.superCluster();
0198 if (sc.isNonnull()) {
0199 eid_sc_E = sc->energy();
0200 eid_sc_eta = sc->eta();
0201 eid_sc_etaWidth = sc->etaWidth();
0202 eid_sc_phiWidth = sc->phiWidth();
0203 sc_Nclus = sc->clustersSize();
0204 }
0205 }
0206
0207
0208 eid_match_seed_dEta = ele.deltaEtaSeedClusterTrackAtCalo();
0209 eid_match_eclu_EoverP = (1. / ele.ecalEnergy()) - (1. / ele.p());
0210 eid_match_SC_EoverP = ele.eSuperClusterOverP();
0211 eid_match_SC_dEta = ele.deltaEtaSuperClusterTrackAtVtx();
0212 eid_match_SC_dPhi = ele.deltaPhiSuperClusterTrackAtVtx();
0213
0214
0215 eid_shape_full5x5_HoverE = ele.full5x5_hcalOverEcal();
0216 eid_shape_full5x5_r9 = ele.full5x5_r9();
0217
0218
0219 eid_rho = rho;
0220
0221 eid_brem_frac = ele.fbrem();
0222 core_shFracHits = ele.shFracInnerHits();
0223
0224
0225 gsf_bdtout1 = unbiased;
0226
0227
0228 if (ele.core().isNonnull()) {
0229 reco::GsfTrackRef gsf = ele.gsfTrack();
0230 if (gsf.isNonnull()) {
0231 reco::SuperClusterRef sc = ele.superCluster();
0232 if (sc.isNonnull()) {
0233
0234 double mass2 = 0.000511 * 0.000511;
0235 float p2 = pow(gsf->p(), 2);
0236 float energy = sqrt(mass2 + p2);
0237 XYZTLorentzVector mom = XYZTLorentzVector(gsf->px(), gsf->py(), gsf->pz(), energy);
0238 XYZTLorentzVector pos = XYZTLorentzVector(gsf->vx(), gsf->vy(), gsf->vz(), 0.);
0239 BaseParticlePropagator propagator(RawParticle(mom, pos, gsf->charge()), 0, 0, field_z);
0240 propagator.propagateToEcalEntrance(true);
0241 bool reach_ECAL = propagator.getSuccess();
0242
0243 GlobalPoint ecal_pos(propagator.particle().x(), propagator.particle().y(), propagator.particle().z());
0244
0245
0246 sc_clus1_nxtal = -999;
0247 sc_clus1_dphi = -999.;
0248 sc_clus2_dphi = -999.;
0249 sc_clus1_deta = -999.;
0250 sc_clus2_deta = -999.;
0251 sc_clus1_E = -999.;
0252 sc_clus2_E = -999.;
0253 sc_clus1_E_ov_p = -999.;
0254 sc_clus2_E_ov_p = -999.;
0255 trackClusterMatching(*sc,
0256 *gsf,
0257 reach_ECAL,
0258 ecal_pos,
0259 sc_clus1_nxtal,
0260 sc_clus1_dphi,
0261 sc_clus2_dphi,
0262 sc_clus1_deta,
0263 sc_clus2_deta,
0264 sc_clus1_E,
0265 sc_clus2_E,
0266 sc_clus1_E_ov_p,
0267 sc_clus2_E_ov_p);
0268 sc_clus1_nxtal = (int)sc_clus1_nxtal;
0269
0270 }
0271 }
0272 }
0273
0274
0275 eid_rho = std::clamp(eid_rho, 0.f, 100.f);
0276 eid_sc_eta = std::clamp(eid_sc_eta, -5.f, 5.f);
0277 eid_shape_full5x5_r9 = std::clamp(eid_shape_full5x5_r9, 0.f, 2.f);
0278 eid_sc_etaWidth = std::clamp(eid_sc_etaWidth, 0.f, 3.14f);
0279 eid_sc_phiWidth = std::clamp(eid_sc_phiWidth, 0.f, 3.14f);
0280 eid_shape_full5x5_HoverE = std::clamp(eid_shape_full5x5_HoverE, 0.f, 50.f);
0281 eid_trk_nhits = std::clamp(eid_trk_nhits, -1.f, 50.f);
0282 eid_trk_chi2red = std::clamp(eid_trk_chi2red, -1.f, 50.f);
0283 eid_gsf_chi2red = std::clamp(eid_gsf_chi2red, -1.f, 100.f);
0284 if (eid_brem_frac < 0.)
0285 eid_brem_frac = -1.;
0286 if (eid_brem_frac > 1.)
0287 eid_brem_frac = 1.;
0288 eid_gsf_nhits = std::clamp(eid_gsf_nhits, -1.f, 50.f);
0289 eid_match_SC_EoverP = std::clamp(eid_match_SC_EoverP, 0.f, 100.f);
0290 eid_match_eclu_EoverP = std::clamp(eid_match_eclu_EoverP, -1.f, 1.f);
0291 eid_match_SC_dEta = std::clamp(eid_match_SC_dEta, -10.f, 10.f);
0292 eid_match_SC_dPhi = std::clamp(eid_match_SC_dPhi, -3.14f, 3.14f);
0293 eid_match_seed_dEta = std::clamp(eid_match_seed_dEta, -10.f, 10.f);
0294 eid_sc_E = std::clamp(eid_sc_E, 0.f, 1000.f);
0295 eid_trk_p = std::clamp(eid_trk_p, -1.f, 1000.f);
0296 gsf_mode_p = std::clamp(gsf_mode_p, 0.f, 1000.f);
0297 core_shFracHits = std::clamp(core_shFracHits, 0.f, 1.f);
0298 gsf_bdtout1 = std::clamp(gsf_bdtout1, -20.f, 20.f);
0299 if (gsf_dr < 0.)
0300 gsf_dr = 5.;
0301 if (gsf_dr > 5.)
0302 gsf_dr = 5.;
0303 if (trk_dr < 0.)
0304 trk_dr = 5.;
0305 if (trk_dr > 5.)
0306 trk_dr = 5.;
0307 sc_Nclus = std::clamp(sc_Nclus, 0.f, 20.f);
0308 sc_clus1_nxtal = std::clamp(sc_clus1_nxtal, 0.f, 100.f);
0309 sc_clus1_dphi = std::clamp(sc_clus1_dphi, -3.14f, 3.14f);
0310 sc_clus2_dphi = std::clamp(sc_clus2_dphi, -3.14f, 3.14f);
0311 sc_clus1_deta = std::clamp(sc_clus1_deta, -5.f, 5.f);
0312 sc_clus2_deta = std::clamp(sc_clus2_deta, -5.f, 5.f);
0313 sc_clus1_E = std::clamp(sc_clus1_E, 0.f, 1000.f);
0314 sc_clus2_E = std::clamp(sc_clus2_E, 0.f, 1000.f);
0315 if (sc_clus1_E_ov_p < 0.)
0316 sc_clus1_E_ov_p = -1.;
0317 if (sc_clus2_E_ov_p < 0.)
0318 sc_clus2_E_ov_p = -1.;
0319
0320
0321 std::vector<float> output = {eid_rho,
0322 eid_sc_eta,
0323 eid_shape_full5x5_r9,
0324 eid_sc_etaWidth,
0325 eid_sc_phiWidth,
0326 eid_shape_full5x5_HoverE,
0327 eid_trk_nhits,
0328 eid_trk_chi2red,
0329 eid_gsf_chi2red,
0330 eid_brem_frac,
0331 eid_gsf_nhits,
0332 eid_match_SC_EoverP,
0333 eid_match_eclu_EoverP,
0334 eid_match_SC_dEta,
0335 eid_match_SC_dPhi,
0336 eid_match_seed_dEta,
0337 eid_sc_E,
0338 eid_trk_p,
0339 gsf_mode_p,
0340 core_shFracHits,
0341 gsf_bdtout1,
0342 gsf_dr,
0343 trk_dr,
0344 sc_Nclus,
0345 sc_clus1_nxtal,
0346 sc_clus1_dphi,
0347 sc_clus2_dphi,
0348 sc_clus1_deta,
0349 sc_clus2_deta,
0350 sc_clus1_E,
0351 sc_clus2_E,
0352 sc_clus1_E_ov_p,
0353 sc_clus2_E_ov_p};
0354 return output;
0355 }
0356
0357
0358
0359 std::vector<float> features_V0(reco::GsfElectron const& ele, float rho, float unbiased) {
0360 float eid_rho = -999.;
0361 float eid_sc_eta = -999.;
0362 float eid_shape_full5x5_r9 = -999.;
0363 float eid_sc_etaWidth = -999.;
0364 float eid_sc_phiWidth = -999.;
0365 float eid_shape_full5x5_HoverE = -999.;
0366 float eid_trk_nhits = -999.;
0367 float eid_trk_chi2red = -999.;
0368 float eid_gsf_chi2red = -999.;
0369 float eid_brem_frac = -999.;
0370 float eid_gsf_nhits = -999.;
0371 float eid_match_SC_EoverP = -999.;
0372 float eid_match_eclu_EoverP = -999.;
0373 float eid_match_SC_dEta = -999.;
0374 float eid_match_SC_dPhi = -999.;
0375 float eid_match_seed_dEta = -999.;
0376 float eid_sc_E = -999.;
0377 float eid_trk_p = -999.;
0378 float gsf_mode_p = -999.;
0379 float core_shFracHits = -999.;
0380 float gsf_bdtout1 = -999.;
0381 float gsf_dr = -999.;
0382 float trk_dr = -999.;
0383 float sc_Nclus = -999.;
0384 float sc_clus1_nxtal = -999.;
0385 float sc_clus1_dphi = -999.;
0386 float sc_clus2_dphi = -999.;
0387 float sc_clus1_deta = -999.;
0388 float sc_clus2_deta = -999.;
0389 float sc_clus1_E = -999.;
0390 float sc_clus2_E = -999.;
0391 float sc_clus1_E_ov_p = -999.;
0392 float sc_clus2_E_ov_p = -999.;
0393
0394
0395 if (ele.core().isNonnull()) {
0396 const auto& trk = ele.closestCtfTrackRef();
0397 if (trk.isNonnull()) {
0398 eid_trk_p = (float)trk->p();
0399 eid_trk_nhits = (float)trk->found();
0400 eid_trk_chi2red = (float)trk->normalizedChi2();
0401 TVector3 trkTV3(0, 0, 0);
0402 trkTV3.SetPtEtaPhi(trk->pt(), trk->eta(), trk->phi());
0403 TVector3 eleTV3(0, 0, 0);
0404 eleTV3.SetPtEtaPhi(ele.pt(), ele.eta(), ele.phi());
0405 trk_dr = eleTV3.DeltaR(trkTV3);
0406 }
0407 }
0408
0409
0410 if (ele.core().isNonnull()) {
0411 const auto& gsf = ele.core()->gsfTrack();
0412 if (gsf.isNonnull()) {
0413 gsf_mode_p = gsf->pMode();
0414 eid_gsf_nhits = (float)gsf->found();
0415 eid_gsf_chi2red = gsf->normalizedChi2();
0416 TVector3 gsfTV3(0, 0, 0);
0417 gsfTV3.SetPtEtaPhi(gsf->ptMode(), gsf->etaMode(), gsf->phiMode());
0418 TVector3 eleTV3(0, 0, 0);
0419 eleTV3.SetPtEtaPhi(ele.pt(), ele.eta(), ele.phi());
0420 gsf_dr = eleTV3.DeltaR(gsfTV3);
0421 }
0422 }
0423
0424
0425 if (ele.core().isNonnull()) {
0426 const auto& sc = ele.core()->superCluster();
0427 if (sc.isNonnull()) {
0428 eid_sc_E = sc->energy();
0429 eid_sc_eta = sc->eta();
0430 eid_sc_etaWidth = sc->etaWidth();
0431 eid_sc_phiWidth = sc->phiWidth();
0432 sc_Nclus = (float)sc->clustersSize();
0433 }
0434 }
0435
0436
0437 eid_match_seed_dEta = ele.deltaEtaSeedClusterTrackAtCalo();
0438 eid_match_eclu_EoverP = (1. / ele.ecalEnergy()) - (1. / ele.p());
0439 eid_match_SC_EoverP = ele.eSuperClusterOverP();
0440 eid_match_SC_dEta = ele.deltaEtaSuperClusterTrackAtVtx();
0441 eid_match_SC_dPhi = ele.deltaPhiSuperClusterTrackAtVtx();
0442
0443
0444 eid_shape_full5x5_HoverE = ele.full5x5_hcalOverEcal();
0445 eid_shape_full5x5_r9 = ele.full5x5_r9();
0446
0447
0448 eid_rho = rho;
0449
0450 eid_brem_frac = ele.fbrem();
0451 core_shFracHits = (float)ele.shFracInnerHits();
0452
0453
0454 gsf_bdtout1 = unbiased;
0455
0456
0457 if (ele.core().isNonnull()) {
0458 const auto& gsf = ele.core()->gsfTrack();
0459 if (gsf.isNonnull()) {
0460 const auto& sc = ele.core()->superCluster();
0461 if (sc.isNonnull()) {
0462
0463 double mass2 = 0.000511 * 0.000511;
0464 float p2 = pow(gsf->p(), 2);
0465 float energy = sqrt(mass2 + p2);
0466 math::XYZTLorentzVector mom = math::XYZTLorentzVector(gsf->px(), gsf->py(), gsf->pz(), energy);
0467 math::XYZTLorentzVector pos = math::XYZTLorentzVector(gsf->vx(), gsf->vy(), gsf->vz(), 0.);
0468 float field_z = 3.8;
0469 BaseParticlePropagator mypart(RawParticle(mom, pos, gsf->charge()), 0, 0, field_z);
0470 mypart.propagateToEcalEntrance(true);
0471 bool reach_ECAL = mypart.getSuccess();
0472
0473
0474 GlobalPoint ecal_pos(
0475 mypart.particle().vertex().x(), mypart.particle().vertex().y(), mypart.particle().vertex().z());
0476
0477
0478 sc_clus1_nxtal = -999.;
0479 sc_clus1_dphi = -999.;
0480 sc_clus2_dphi = -999.;
0481 sc_clus1_deta = -999.;
0482 sc_clus2_deta = -999.;
0483 sc_clus1_E = -999.;
0484 sc_clus2_E = -999.;
0485 sc_clus1_E_ov_p = -999.;
0486 sc_clus2_E_ov_p = -999.;
0487 trackClusterMatching(*sc,
0488 *gsf,
0489 reach_ECAL,
0490 ecal_pos,
0491 sc_clus1_nxtal,
0492 sc_clus1_dphi,
0493 sc_clus2_dphi,
0494 sc_clus1_deta,
0495 sc_clus2_deta,
0496 sc_clus1_E,
0497 sc_clus2_E,
0498 sc_clus1_E_ov_p,
0499 sc_clus2_E_ov_p);
0500
0501 }
0502 }
0503 }
0504
0505
0506 eid_sc_eta = std::clamp(eid_sc_eta, -5.f, 5.f);
0507 eid_shape_full5x5_r9 = std::clamp(eid_shape_full5x5_r9, 0.f, 2.f);
0508 eid_sc_etaWidth = std::clamp(eid_sc_etaWidth, 0.f, 3.14f);
0509 eid_sc_phiWidth = std::clamp(eid_sc_phiWidth, 0.f, 3.14f);
0510 eid_shape_full5x5_HoverE = std::clamp(eid_shape_full5x5_HoverE, 0.f, 50.f);
0511 eid_trk_nhits = std::clamp(eid_trk_nhits, -1.f, 50.f);
0512 eid_trk_chi2red = std::clamp(eid_trk_chi2red, -1.f, 50.f);
0513 eid_gsf_chi2red = std::clamp(eid_gsf_chi2red, -1.f, 100.f);
0514 if (eid_brem_frac < 0.)
0515 eid_brem_frac = -1.;
0516 if (eid_brem_frac > 1.)
0517 eid_brem_frac = 1.;
0518 eid_gsf_nhits = std::clamp(eid_gsf_nhits, -1.f, 50.f);
0519 eid_match_SC_EoverP = std::clamp(eid_match_SC_EoverP, 0.f, 100.f);
0520 eid_match_eclu_EoverP = std::clamp(eid_match_eclu_EoverP, -1.f, 1.f);
0521 eid_match_SC_dEta = std::clamp(eid_match_SC_dEta, -10.f, 10.f);
0522 eid_match_SC_dPhi = std::clamp(eid_match_SC_dPhi, -3.14f, 3.14f);
0523 eid_match_seed_dEta = std::clamp(eid_match_seed_dEta, -10.f, 10.f);
0524 eid_sc_E = std::clamp(eid_sc_E, 0.f, 1000.f);
0525 eid_trk_p = std::clamp(eid_trk_p, -1.f, 1000.f);
0526 gsf_mode_p = std::clamp(gsf_mode_p, 0.f, 1000.f);
0527 core_shFracHits = std::clamp(core_shFracHits, 0.f, 1.f);
0528 gsf_bdtout1 = std::clamp(gsf_bdtout1, -20.f, 20.f);
0529 if (gsf_dr < 0.)
0530 gsf_dr = 5.;
0531 if (gsf_dr > 5.)
0532 gsf_dr = 5.;
0533 if (trk_dr < 0.)
0534 trk_dr = 5.;
0535 if (trk_dr > 5.)
0536 trk_dr = 5.;
0537 sc_Nclus = std::clamp(sc_Nclus, 0.f, 20.f);
0538 sc_clus1_nxtal = std::clamp(sc_clus1_nxtal, 0.f, 100.f);
0539 if (sc_clus1_dphi < -3.14)
0540 sc_clus1_dphi = -5.;
0541 if (sc_clus1_dphi > 3.14)
0542 sc_clus1_dphi = 5.;
0543 if (sc_clus2_dphi < -3.14)
0544 sc_clus2_dphi = -5.;
0545 if (sc_clus2_dphi > 3.14)
0546 sc_clus2_dphi = 5.;
0547 sc_clus1_deta = std::clamp(sc_clus1_deta, -5.f, 5.f);
0548 sc_clus2_deta = std::clamp(sc_clus2_deta, -5.f, 5.f);
0549 sc_clus1_E = std::clamp(sc_clus1_E, 0.f, 1000.f);
0550 sc_clus2_E = std::clamp(sc_clus2_E, 0.f, 1000.f);
0551 if (sc_clus1_E_ov_p < 0.)
0552 sc_clus1_E_ov_p = -1.;
0553 if (sc_clus2_E_ov_p < 0.)
0554 sc_clus2_E_ov_p = -1.;
0555
0556
0557 std::vector<float> output = {eid_rho,
0558 eid_sc_eta,
0559 eid_shape_full5x5_r9,
0560 eid_sc_etaWidth,
0561 eid_sc_phiWidth,
0562 eid_shape_full5x5_HoverE,
0563 eid_trk_nhits,
0564 eid_trk_chi2red,
0565 eid_gsf_chi2red,
0566 eid_brem_frac,
0567 eid_gsf_nhits,
0568 eid_match_SC_EoverP,
0569 eid_match_eclu_EoverP,
0570 eid_match_SC_dEta,
0571 eid_match_SC_dPhi,
0572 eid_match_seed_dEta,
0573 eid_sc_E,
0574 eid_trk_p,
0575 gsf_mode_p,
0576 core_shFracHits,
0577 gsf_bdtout1,
0578 gsf_dr,
0579 trk_dr,
0580 sc_Nclus,
0581 sc_clus1_nxtal,
0582 sc_clus1_dphi,
0583 sc_clus2_dphi,
0584 sc_clus1_deta,
0585 sc_clus2_deta,
0586 sc_clus1_E,
0587 sc_clus2_E,
0588 sc_clus1_E_ov_p,
0589 sc_clus2_E_ov_p};
0590 return output;
0591 }
0592
0593
0594
0595 void findEnergeticClusters(
0596 reco::SuperCluster const& sc, int& clusNum, float& maxEne1, float& maxEne2, int& i1, int& i2) {
0597 if (sc.clustersSize() > 0 && sc.clustersBegin() != sc.clustersEnd()) {
0598 for (auto const& cluster : sc.clusters()) {
0599 if (cluster->energy() > maxEne1) {
0600 maxEne1 = cluster->energy();
0601 i1 = clusNum;
0602 }
0603 clusNum++;
0604 }
0605 if (sc.clustersSize() > 1) {
0606 clusNum = 0;
0607 for (auto const& cluster : sc.clusters()) {
0608 if (clusNum != i1) {
0609 if (cluster->energy() > maxEne2) {
0610 maxEne2 = cluster->energy();
0611 i2 = clusNum;
0612 }
0613 }
0614 clusNum++;
0615 }
0616 }
0617 }
0618 }
0619
0620
0621
0622 void trackClusterMatching(reco::SuperCluster const& sc,
0623 reco::GsfTrack const& gsf,
0624 bool const& reach_ECAL,
0625 GlobalPoint const& ecal_pos,
0626 float& sc_clus1_nxtal,
0627 float& sc_clus1_dphi,
0628 float& sc_clus2_dphi,
0629 float& sc_clus1_deta,
0630 float& sc_clus2_deta,
0631 float& sc_clus1_E,
0632 float& sc_clus2_E,
0633 float& sc_clus1_E_ov_p,
0634 float& sc_clus2_E_ov_p) {
0635
0636 int clusNum = 0;
0637 float maxEne1 = -1;
0638 float maxEne2 = -1;
0639 int i1 = -1;
0640 int i2 = -1;
0641 findEnergeticClusters(sc, clusNum, maxEne1, maxEne2, i1, i2);
0642
0643
0644 clusNum = 0;
0645 if (sc.clustersSize() > 0 && sc.clustersBegin() != sc.clustersEnd()) {
0646 for (auto const& cluster : sc.clusters()) {
0647 float deta = ecal_pos.eta() - cluster->eta();
0648 float dphi = reco::deltaPhi(ecal_pos.phi(), cluster->phi());
0649 if (clusNum == i1) {
0650 sc_clus1_E = cluster->energy();
0651 if (gsf.pMode() > 0)
0652 sc_clus1_E_ov_p = cluster->energy() / gsf.pMode();
0653 sc_clus1_nxtal = (float)cluster->size();
0654 if (reach_ECAL > 0) {
0655 sc_clus1_deta = deta;
0656 sc_clus1_dphi = dphi;
0657 }
0658 } else if (clusNum == i2) {
0659 sc_clus2_E = cluster->energy();
0660 if (gsf.pMode() > 0)
0661 sc_clus2_E_ov_p = cluster->energy() / gsf.pMode();
0662 if (reach_ECAL > 0) {
0663 sc_clus2_deta = deta;
0664 sc_clus2_dphi = dphi;
0665 }
0666 }
0667 clusNum++;
0668 }
0669 }
0670 }
0671
0672 }
