Back to home page

Project CMSSW displayed by LXR

 
 

    


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.;  // must be last (not used by unbiased model)
0047 
0048     // Tracks
0049     const auto& trk = ecal.trackRef();  // reco::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);  //@@ to be consistent with the training based on 94X MC
0060       }
0061       ktf_ecal_cluster_dphi_ *= trk->charge();  //@@ to be consistent with the training based on 94X MC
0062     }
0063 
0064     // Rho
0065     rho_ = static_cast<float>(rho);
0066 
0067     // ECAL clusters
0068     const auto& ecal_clu = ecal.clusterRef();  // reco::PFClusterRef
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     // HCAL clusters
0090     const auto& hcal_clu = hcal.clusterRef();  // reco::PFClusterRef
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     // PreId
0098     preid_gsf_dpt_ = ecal.dpt();
0099     preid_trk_gsf_chiratio_ = ecal.chi2Ratio();
0100     preid_gsf_chi2red_ = ecal.gsfChi2();
0101 
0102     // Set contents of vector
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 }  // namespace lowptgsfeleseed
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     // KF tracks
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     // GSF tracks
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     // Super clusters
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     // Track-cluster matching
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     // Shower shape vars
0215     eid_shape_full5x5_HoverE = ele.full5x5_hcalOverEcal();
0216     eid_shape_full5x5_r9 = ele.full5x5_r9();
0217 
0218     // Misc
0219     eid_rho = rho;
0220 
0221     eid_brem_frac = ele.fbrem();
0222     core_shFracHits = ele.shFracInnerHits();
0223 
0224     // Unbiased BDT from ElectronSeed
0225     gsf_bdtout1 = unbiased;
0226 
0227     // Clusters
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           // Propagate electron track to ECAL surface
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);   // true only first half loop , false more than one loop
0241           bool reach_ECAL = propagator.getSuccess();  // 0 does not reach ECAL, 1 yes barrel, 2 yes endcaps
0242                                                       // ECAL entry point for track
0243           GlobalPoint ecal_pos(propagator.particle().x(), propagator.particle().y(), propagator.particle().z());
0244 
0245           // Track-cluster matching for most energetic clusters
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         }  // sc.isNonnull()
0271       }    // gsf.isNonnull()
0272     }      // clusters
0273 
0274     // Out-of-range
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     // Set contents of vector
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   // feature list for original models (2019Aug07 and earlier)
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     // KF tracks
0395     if (ele.core().isNonnull()) {
0396       const auto& trk = ele.closestCtfTrackRef();  // reco::TrackRef
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     // GSF tracks
0410     if (ele.core().isNonnull()) {
0411       const auto& gsf = ele.core()->gsfTrack();  // reco::GsfTrackRef
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     // Super clusters
0425     if (ele.core().isNonnull()) {
0426       const auto& sc = ele.core()->superCluster();  // reco::SuperClusterRef
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     // Track-cluster matching
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     // Shower shape vars
0444     eid_shape_full5x5_HoverE = ele.full5x5_hcalOverEcal();
0445     eid_shape_full5x5_r9 = ele.full5x5_r9();
0446 
0447     // Misc
0448     eid_rho = rho;
0449 
0450     eid_brem_frac = ele.fbrem();
0451     core_shFracHits = (float)ele.shFracInnerHits();
0452 
0453     // Unbiased BDT from ElectronSeed
0454     gsf_bdtout1 = unbiased;
0455 
0456     // Clusters
0457     if (ele.core().isNonnull()) {
0458       const auto& gsf = ele.core()->gsfTrack();  // reco::GsfTrackRef
0459       if (gsf.isNonnull()) {
0460         const auto& sc = ele.core()->superCluster();  // reco::SuperClusterRef
0461         if (sc.isNonnull()) {
0462           // Propagate electron track to ECAL surface
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);   // true only first half loop , false more than one loop
0471           bool reach_ECAL = mypart.getSuccess();  // 0 does not reach ECAL, 1 yes barrel, 2 yes endcaps
0472 
0473           // ECAL entry point for track
0474           GlobalPoint ecal_pos(
0475               mypart.particle().vertex().x(), mypart.particle().vertex().y(), mypart.particle().vertex().z());
0476 
0477           // Track-cluster matching for most energetic clusters
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         }  // sc.isNonnull()
0502       }    // gsf.isNonnull()
0503     }      // clusters
0504 
0505     // Out-of-range
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     // Set contents of vector
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   // Find most energetic clusters
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     }  // loop over clusters
0618   }
0619 
0620   ////////////////////////////////////////////////////////////////////////////////
0621   // Track-cluster matching for most energetic clusters
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     // Iterate through ECAL clusters and sort in energy
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     // track-clusters match
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 }  // namespace lowptgsfeleid