Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-11-05 00:08:11

0001 #include "RecoParticleFlow/PFProducer/interface/PFEGammaAlgo.h"
0002 #include "RecoParticleFlow/PFProducer/interface/PFMuonAlgo.h"
0003 #include "DataFormats/ParticleFlowReco/interface/PFRecHitFraction.h"
0004 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
0005 #include "RecoParticleFlow/PFClusterTools/interface/ClusterClusterMapping.h"
0006 #include "DataFormats/ParticleFlowReco/interface/PFLayer.h"
0007 #include "DataFormats/ParticleFlowReco/interface/PFRecHit.h"
0008 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0009 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0010 #include "DataFormats/EcalDetId/interface/ESDetId.h"
0011 #include "DataFormats/EgammaCandidates/interface/Photon.h"
0012 #include "RecoParticleFlow/PFClusterTools/interface/PFEnergyCalibration.h"
0013 #include "RecoParticleFlow/PFClusterTools/interface/PFEnergyResolution.h"
0014 #include "CommonTools/ParticleFlow/interface/PFClusterWidthAlgo.h"
0015 #include "RecoParticleFlow/PFTracking/interface/PFTrackAlgoTools.h"
0016 #include "DataFormats/Common/interface/RefToPtr.h"
0017 #include "RecoEcal/EgammaCoreTools/interface/Mustache.h"
0018 #include "DataFormats/Math/interface/deltaPhi.h"
0019 #include "DataFormats/Math/interface/deltaR.h"
0020 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0021 
0022 #include <TFile.h>
0023 #include <TVector2.h>
0024 #include <iomanip>
0025 #include <algorithm>
0026 #include <functional>
0027 #include <numeric>
0028 #include <TMath.h>
0029 
0030 // include combinations header (was never incorporated in boost)
0031 #include "combination.hpp"
0032 
0033 // just for now do this
0034 //#define PFLOW_DEBUG
0035 
0036 #ifdef PFLOW_DEBUG
0037 #define docast(x, y) dynamic_cast<x>(y)
0038 #define LOGVERB(x) edm::LogVerbatim(x)
0039 #define LOGWARN(x) edm::LogWarning(x)
0040 #define LOGERR(x) edm::LogError(x)
0041 #define LOGDRESSED(x) edm::LogInfo(x)
0042 #else
0043 #define docast(x, y) static_cast<x>(y)
0044 #define LOGVERB(x) LogTrace(x)
0045 #define LOGWARN(x) edm::LogWarning(x)
0046 #define LOGERR(x) edm::LogError(x)
0047 #define LOGDRESSED(x) LogDebug(x)
0048 #endif
0049 
0050 using namespace std;
0051 using namespace reco;
0052 using namespace std::placeholders;  // for _1, _2, _3...
0053 
0054 namespace {
0055   typedef PFEGammaAlgo::PFSCElement SCElement;
0056   typedef PFEGammaAlgo::EEtoPSAssociation EEtoPSAssociation;
0057   typedef std::pair<CaloClusterPtr::key_type, CaloClusterPtr> EEtoPSElement;
0058   typedef PFEGammaAlgo::PFClusterElement ClusterElement;
0059 
0060   class SeedMatchesToProtoObject {
0061   public:
0062     SeedMatchesToProtoObject(const reco::ElectronSeedRef& s)
0063         : scfromseed_(s->caloCluster().castTo<reco::SuperClusterRef>()) {
0064       ispfsc_ = false;
0065       if (scfromseed_.isNonnull()) {
0066         const edm::Ptr<reco::PFCluster> testCast(scfromseed_->seed());
0067         ispfsc_ = testCast.isNonnull();
0068       }
0069     }
0070     bool operator()(const PFEGammaAlgo::ProtoEGObject& po) {
0071       if (scfromseed_.isNull() || !po.parentSC)
0072         return false;
0073       if (ispfsc_) {
0074         return (scfromseed_->seed() == po.parentSC->superClusterRef()->seed());
0075       }
0076       return (scfromseed_->seed()->seed() == po.parentSC->superClusterRef()->seed()->seed());
0077     }
0078 
0079   private:
0080     const reco::SuperClusterRef scfromseed_;
0081     bool ispfsc_;
0082   };
0083 
0084   template <bool useConvs = false>
0085   bool elementNotCloserToOther(const reco::PFBlockRef& block,
0086                                const PFBlockElement::Type& keytype,
0087                                const size_t key,
0088                                const PFBlockElement::Type& valtype,
0089                                const size_t test,
0090                                const float EoPin_cut = 1.0e6) {
0091     constexpr reco::PFBlockElement::TrackType ConvType = reco::PFBlockElement::T_FROM_GAMMACONV;
0092     auto elemkey = &(block->elements()[key]);
0093     // this is inside out but I just want something that works right now
0094     switch (elemkey->type()) {
0095       case reco::PFBlockElement::GSF: {
0096         const reco::PFBlockElementGsfTrack* elemasgsf = docast(const reco::PFBlockElementGsfTrack*, elemkey);
0097         if (elemasgsf && valtype == PFBlockElement::ECAL) {
0098           const ClusterElement* elemasclus = static_cast<const ClusterElement*>(&(block->elements()[test]));
0099           float cluster_e = elemasclus->clusterRef()->correctedEnergy();
0100           float trk_pin = elemasgsf->Pin().P();
0101           if (cluster_e / trk_pin > EoPin_cut) {
0102             LOGDRESSED("elementNotCloserToOther") << "GSF track failed EoP cut to match with cluster!";
0103             return false;
0104           }
0105         }
0106       } break;
0107       case reco::PFBlockElement::TRACK: {
0108         const reco::PFBlockElementTrack* elemaskf = docast(const reco::PFBlockElementTrack*, elemkey);
0109         if (elemaskf && valtype == PFBlockElement::ECAL) {
0110           const ClusterElement* elemasclus = static_cast<const ClusterElement*>(&(block->elements()[test]));
0111           float cluster_e = elemasclus->clusterRef()->correctedEnergy();
0112           float trk_pin = std::sqrt(elemaskf->trackRef()->innerMomentum().mag2());
0113           if (cluster_e / trk_pin > EoPin_cut) {
0114             LOGDRESSED("elementNotCloserToOther") << "KF track failed EoP cut to match with cluster!";
0115             return false;
0116           }
0117         }
0118       } break;
0119       default:
0120         break;
0121     }
0122 
0123     const float dist = block->dist(key, test, block->linkData(), reco::PFBlock::LINKTEST_ALL);
0124     if (dist == -1.0f)
0125       return false;  // don't associate non-linked elems
0126     std::multimap<double, unsigned> dists_to_val;
0127     block->associatedElements(test, block->linkData(), dists_to_val, keytype, reco::PFBlock::LINKTEST_ALL);
0128 
0129     for (const auto& valdist : dists_to_val) {
0130       const size_t idx = valdist.second;
0131       // check track types for conversion info
0132       switch (keytype) {
0133         case reco::PFBlockElement::GSF: {
0134           const reco::PFBlockElementGsfTrack* elemasgsf =
0135               docast(const reco::PFBlockElementGsfTrack*, &(block->elements()[idx]));
0136           if (!useConvs && elemasgsf->trackType(ConvType))
0137             return false;
0138           if (elemasgsf && valtype == PFBlockElement::ECAL) {
0139             const ClusterElement* elemasclus = docast(const ClusterElement*, &(block->elements()[test]));
0140             float cluster_e = elemasclus->clusterRef()->correctedEnergy();
0141             float trk_pin = elemasgsf->Pin().P();
0142             if (cluster_e / trk_pin > EoPin_cut)
0143               continue;
0144           }
0145         } break;
0146         case reco::PFBlockElement::TRACK: {
0147           const reco::PFBlockElementTrack* elemaskf =
0148               docast(const reco::PFBlockElementTrack*, &(block->elements()[idx]));
0149           if (!useConvs && elemaskf->trackType(ConvType))
0150             return false;
0151           if (elemaskf && valtype == PFBlockElement::ECAL) {
0152             const ClusterElement* elemasclus = static_cast<const ClusterElement*>(&(block->elements()[test]));
0153             float cluster_e = elemasclus->clusterRef()->correctedEnergy();
0154             float trk_pin = std::sqrt(elemaskf->trackRef()->innerMomentum().mag2());
0155             if (cluster_e / trk_pin > EoPin_cut)
0156               continue;
0157           }
0158         } break;
0159         default:
0160           break;
0161       }
0162       if (valdist.first < dist && idx != key) {
0163         LOGDRESSED("elementNotCloserToOther")
0164             << "key element of type " << keytype << " is closer to another element of type" << valtype << std::endl;
0165         return false;  // false if closer element of specified type found
0166       }
0167     }
0168     return true;
0169   }
0170 
0171   template <class Element1, class Element2>
0172   bool compatibleEoPOut(const Element1& e, const Element2& comp) {
0173     if (PFBlockElement::ECAL != e.type()) {
0174       return false;
0175     }
0176     const ClusterElement& elemascluster = docast(ClusterElement const&, e);
0177     const float gsf_eta_diff = std::abs(comp.positionAtECALEntrance().eta() - comp.Pout().eta());
0178     const reco::PFClusterRef& cRef = elemascluster.clusterRef();
0179     return (gsf_eta_diff <= 0.3 && cRef->energy() / comp.Pout().t() <= 5);
0180   }
0181 
0182   constexpr reco::PFBlockElement::TrackType ConvType = reco::PFBlockElement::T_FROM_GAMMACONV;
0183 
0184   template <PFBlockElement::Type keytype, PFBlockElement::Type valtype, bool useConv = false>
0185 
0186   struct NotCloserToOther {
0187     const reco::PFBlockElement* comp;
0188     const reco::PFBlockRef& block;
0189     const float EoPin_cut;
0190     NotCloserToOther(const reco::PFBlockRef& b, const reco::PFBlockElement* e, const float EoPcut = 1.0e6)
0191         : comp(e), block(b), EoPin_cut(EoPcut) {}
0192     template <class T>
0193     bool operator()(const T& e) {
0194       if (!e.flag() || valtype != e->type())
0195         return false;
0196       return elementNotCloserToOther<useConv>(block, keytype, comp->index(), valtype, e->index(), EoPin_cut);
0197     }
0198   };
0199 
0200   struct LesserByDistance {
0201     const reco::PFBlockElement* comp;
0202     const reco::PFBlockRef& block;
0203     const reco::PFBlock::LinkData& links;
0204     LesserByDistance(const reco::PFBlockRef& b, const reco::PFBlock::LinkData& l, const reco::PFBlockElement* e)
0205         : comp(e), block(b), links(l) {}
0206     bool operator()(FlaggedPtr<const reco::PFBlockElement> const& e1,
0207                     FlaggedPtr<const reco::PFBlockElement> const& e2) {
0208       double dist1 = block->dist(comp->index(), e1->index(), links, reco::PFBlock::LINKTEST_ALL);
0209       double dist2 = block->dist(comp->index(), e2->index(), links, reco::PFBlock::LINKTEST_ALL);
0210       dist1 = (dist1 == -1.0 ? 1e6 : dist1);
0211       dist2 = (dist2 == -1.0 ? 1e6 : dist2);
0212       return dist1 < dist2;
0213     }
0214   };
0215 
0216   bool isROLinkedByClusterOrTrack(const PFEGammaAlgo::ProtoEGObject& RO1, const PFEGammaAlgo::ProtoEGObject& RO2) {
0217     // also don't allow ROs where both have clusters
0218     // and GSF tracks to merge (10 Dec 2013)
0219     if (!RO1.primaryGSFs.empty() && !RO2.primaryGSFs.empty()) {
0220       LOGDRESSED("isROLinkedByClusterOrTrack") << "cannot merge, both have GSFs!" << std::endl;
0221       return false;
0222     }
0223     // don't allow EB/EE to mix (11 Sept 2013)
0224     if (!RO1.ecalclusters.empty() && !RO2.ecalclusters.empty()) {
0225       if (RO1.ecalclusters.front()->clusterRef()->layer() != RO2.ecalclusters.front()->clusterRef()->layer()) {
0226         LOGDRESSED("isROLinkedByClusterOrTrack") << "cannot merge, different ECAL types!" << std::endl;
0227         return false;
0228       }
0229     }
0230     const reco::PFBlockRef& blk = RO1.parentBlock;
0231     bool not_closer;
0232     // check links track -> cluster
0233     for (const auto& cluster : RO1.ecalclusters) {
0234       for (const auto& primgsf : RO2.primaryGSFs) {
0235         not_closer = elementNotCloserToOther(blk, cluster->type(), cluster->index(), primgsf->type(), primgsf->index());
0236         if (not_closer) {
0237           LOGDRESSED("isROLinkedByClusterOrTrack") << "merged by cluster to primary GSF" << std::endl;
0238           return true;
0239         } else {
0240           LOGDRESSED("isROLinkedByClusterOrTrack") << "cluster to primary GSF failed since"
0241                                                    << " cluster closer to another GSF" << std::endl;
0242         }
0243       }
0244       for (const auto& primkf : RO2.primaryKFs) {
0245         not_closer = elementNotCloserToOther(blk, cluster->type(), cluster->index(), primkf->type(), primkf->index());
0246         if (not_closer) {
0247           LOGDRESSED("isROLinkedByClusterOrTrack") << "merged by cluster to primary KF" << std::endl;
0248           return true;
0249         }
0250       }
0251       for (const auto& secdkf : RO2.secondaryKFs) {
0252         not_closer = elementNotCloserToOther(blk, cluster->type(), cluster->index(), secdkf->type(), secdkf->index());
0253         if (not_closer) {
0254           LOGDRESSED("isROLinkedByClusterOrTrack") << "merged by cluster to secondary KF" << std::endl;
0255           return true;
0256         }
0257       }
0258       // check links brem -> cluster
0259       for (const auto& brem : RO2.brems) {
0260         not_closer = elementNotCloserToOther(blk, cluster->type(), cluster->index(), brem->type(), brem->index());
0261         if (not_closer) {
0262           LOGDRESSED("isROLinkedByClusterOrTrack") << "merged by cluster to brem KF" << std::endl;
0263           return true;
0264         }
0265       }
0266     }
0267     // check links primary gsf -> secondary kf
0268     for (const auto& primgsf : RO1.primaryGSFs) {
0269       for (const auto& secdkf : RO2.secondaryKFs) {
0270         not_closer = elementNotCloserToOther(blk, primgsf->type(), primgsf->index(), secdkf->type(), secdkf->index());
0271         if (not_closer) {
0272           LOGDRESSED("isROLinkedByClusterOrTrack") << "merged by GSF to secondary KF" << std::endl;
0273           return true;
0274         }
0275       }
0276     }
0277     // check links primary kf -> secondary kf
0278     for (const auto& primkf : RO1.primaryKFs) {
0279       for (const auto& secdkf : RO2.secondaryKFs) {
0280         not_closer = elementNotCloserToOther(blk, primkf->type(), primkf->index(), secdkf->type(), secdkf->index());
0281         if (not_closer) {
0282           LOGDRESSED("isROLinkedByClusterOrTrack") << "merged by primary KF to secondary KF" << std::endl;
0283           return true;
0284         }
0285       }
0286     }
0287     // check links secondary kf -> secondary kf
0288     for (const auto& secdkf1 : RO1.secondaryKFs) {
0289       for (const auto& secdkf2 : RO2.secondaryKFs) {
0290         not_closer =
0291             elementNotCloserToOther<true>(blk, secdkf1->type(), secdkf1->index(), secdkf2->type(), secdkf2->index());
0292         if (not_closer) {
0293           LOGDRESSED("isROLinkedByClusterOrTrack") << "merged by secondary KF to secondary KF" << std::endl;
0294           return true;
0295         }
0296       }
0297     }
0298     return false;
0299   }
0300 
0301   bool testIfROMergableByLink(const PFEGammaAlgo::ProtoEGObject& ro, PFEGammaAlgo::ProtoEGObject& comp) {
0302     const bool result = (isROLinkedByClusterOrTrack(comp, ro) || isROLinkedByClusterOrTrack(ro, comp));
0303     return result;
0304   }
0305 
0306   std::vector<const ClusterElement*> getSCAssociatedECALsSafe(
0307       const reco::SuperClusterRef& scref, std::vector<FlaggedPtr<const reco::PFBlockElement>>& ecals) {
0308     std::vector<const ClusterElement*> cluster_list;
0309     auto sccl = scref->clustersBegin();
0310     auto scend = scref->clustersEnd();
0311     auto pfc = ecals.begin();
0312     auto pfcend = ecals.end();
0313     for (; sccl != scend; ++sccl) {
0314       std::vector<const ClusterElement*> matched_pfcs;
0315       const double eg_energy = (*sccl)->energy();
0316 
0317       for (pfc = ecals.begin(); pfc != pfcend; ++pfc) {
0318         const ClusterElement* pfcel = docast(const ClusterElement*, pfc->get());
0319         const bool matched = ClusterClusterMapping::overlap(**sccl, *(pfcel->clusterRef()));
0320         // need to protect against high energy clusters being attached
0321         // to low-energy SCs
0322         if (matched && pfcel->clusterRef()->energy() < 1.2 * scref->energy()) {
0323           matched_pfcs.push_back(pfcel);
0324         }
0325       }
0326       std::sort(matched_pfcs.begin(), matched_pfcs.end());
0327 
0328       double min_residual = 1e6;
0329       std::vector<const ClusterElement*> best_comb;
0330       for (size_t i = 1; i <= matched_pfcs.size(); ++i) {
0331         //now we find the combination of PF-clusters which
0332         //has the smallest energy residual with respect to the
0333         //EG-cluster we are looking at now
0334         do {
0335           double energy = std::accumulate(
0336               matched_pfcs.begin(), matched_pfcs.begin() + i, 0.0, [](const double a, const ClusterElement* c) {
0337                 return a + c->clusterRef()->energy();
0338               });
0339           const double resid = std::abs(energy - eg_energy);
0340           if (resid < min_residual) {
0341             best_comb.clear();
0342             best_comb.reserve(i);
0343             min_residual = resid;
0344             best_comb.insert(best_comb.begin(), matched_pfcs.begin(), matched_pfcs.begin() + i);
0345           }
0346         } while (notboost::next_combination(matched_pfcs.begin(), matched_pfcs.begin() + i, matched_pfcs.end()));
0347       }
0348       for (const auto& clelem : best_comb) {
0349         if (std::find(cluster_list.begin(), cluster_list.end(), clelem) == cluster_list.end()) {
0350           cluster_list.push_back(clelem);
0351         }
0352       }
0353     }
0354     return cluster_list;
0355   }
0356   bool addPFClusterToROSafe(const ClusterElement* cl, PFEGammaAlgo::ProtoEGObject& RO) {
0357     if (RO.ecalclusters.empty()) {
0358       RO.ecalclusters.emplace_back(cl, true);
0359       return true;
0360     } else {
0361       const PFLayer::Layer clayer = cl->clusterRef()->layer();
0362       const PFLayer::Layer blayer = RO.ecalclusters.back()->clusterRef()->layer();
0363       if (clayer == blayer) {
0364         RO.ecalclusters.emplace_back(cl, true);
0365         return true;
0366       }
0367     }
0368     return false;
0369   }
0370 
0371   // sets the cluster best associated to the GSF track
0372   // leave it null if no GSF track
0373   void setROElectronCluster(PFEGammaAlgo::ProtoEGObject& RO) {
0374     if (RO.ecalclusters.empty())
0375       return;
0376     RO.lateBrem = -1;
0377     RO.firstBrem = -1;
0378     RO.nBremsWithClusters = -1;
0379     const reco::PFBlockElementBrem *firstBrem = nullptr, *lastBrem = nullptr;
0380     const reco::PFBlockElementCluster *bremCluster = nullptr, *gsfCluster = nullptr, *kfCluster = nullptr,
0381                                       *gsfCluster_noassc = nullptr;
0382     const reco::PFBlockRef& parent = RO.parentBlock;
0383     int nBremClusters = 0;
0384     constexpr float maxDist = 1e6;
0385     float mDist_gsf(maxDist), mDist_gsf_noassc(maxDist), mDist_kf(maxDist);
0386     for (const auto& cluster : RO.ecalclusters) {
0387       for (const auto& gsf : RO.primaryGSFs) {
0388         const bool hasclu =
0389             elementNotCloserToOther(parent, gsf->type(), gsf->index(), cluster->type(), cluster->index());
0390         const float deta = std::abs(cluster->clusterRef()->positionREP().eta() - gsf->positionAtECALEntrance().eta());
0391         const float dphi = std::abs(
0392             TVector2::Phi_mpi_pi(cluster->clusterRef()->positionREP().phi() - gsf->positionAtECALEntrance().phi()));
0393         const float dist = std::hypot(deta, dphi);
0394         if (hasclu && dist < mDist_gsf) {
0395           gsfCluster = cluster.get();
0396           mDist_gsf = dist;
0397         } else if (dist < mDist_gsf_noassc) {
0398           gsfCluster_noassc = cluster.get();
0399           mDist_gsf_noassc = dist;
0400         }
0401       }
0402       for (const auto& kf : RO.primaryKFs) {
0403         const bool hasclu = elementNotCloserToOther(parent, kf->type(), kf->index(), cluster->type(), cluster->index());
0404         const float dist = parent->dist(cluster->index(), kf->index(), parent->linkData(), reco::PFBlock::LINKTEST_ALL);
0405         if (hasclu && dist < mDist_kf) {
0406           kfCluster = cluster.get();
0407           mDist_kf = dist;
0408         }
0409       }
0410       for (const auto& brem : RO.brems) {
0411         const bool hasclu =
0412             elementNotCloserToOther(parent, brem->type(), brem->index(), cluster->type(), cluster->index());
0413         if (hasclu) {
0414           ++nBremClusters;
0415           if (!firstBrem || (firstBrem->indTrajPoint() - 2 > brem->indTrajPoint() - 2)) {
0416             firstBrem = brem;
0417           }
0418           if (!lastBrem || (lastBrem->indTrajPoint() - 2 < brem->indTrajPoint() - 2)) {
0419             lastBrem = brem;
0420             bremCluster = cluster.get();
0421           }
0422         }
0423       }
0424     }
0425     if (!gsfCluster && !kfCluster && !bremCluster) {
0426       gsfCluster = gsfCluster_noassc;
0427     }
0428     RO.nBremsWithClusters = nBremClusters;
0429     RO.lateBrem = 0;
0430     if (gsfCluster) {
0431       RO.electronClusters.push_back(gsfCluster);
0432     } else if (kfCluster) {
0433       RO.electronClusters.push_back(kfCluster);
0434     }
0435     if (bremCluster && !gsfCluster && !kfCluster) {
0436       RO.electronClusters.push_back(bremCluster);
0437     }
0438     if (firstBrem && RO.ecalclusters.size() > 1) {
0439       RO.firstBrem = firstBrem->indTrajPoint() - 2;
0440       if (bremCluster == gsfCluster)
0441         RO.lateBrem = 1;
0442     }
0443   }
0444 }  // namespace
0445 
0446 PFEGammaAlgo::PFEGammaAlgo(const PFEGammaAlgo::PFEGConfigInfo& cfg,
0447                            GBRForests const& gbrForests,
0448                            EEtoPSAssociation const& eetops,
0449                            ESEEIntercalibConstants const& esEEInterCalib,
0450                            ESChannelStatus const& channelStatus,
0451                            reco::Vertex const& primaryVertex)
0452     : gbrForests_(gbrForests),
0453       eetops_(eetops),
0454       cfg_(cfg),
0455       primaryVertex_(primaryVertex),
0456       channelStatus_(channelStatus) {
0457   thePFEnergyCalibration_.initAlphaGamma_ESplanes_fromDB(&esEEInterCalib);
0458 }
0459 
0460 float PFEGammaAlgo::evaluateSingleLegMVA(const reco::PFBlockRef& blockRef,
0461                                          const reco::Vertex& primaryVtx,
0462                                          unsigned int trackIndex) {
0463   const reco::PFBlock& block = *blockRef;
0464   const edm::OwnVector<reco::PFBlockElement>& elements = block.elements();
0465   //use this to store linkdata in the associatedElements function below
0466   const PFBlock::LinkData& linkData = block.linkData();
0467   //calculate MVA Variables
0468   const float chi2 = elements[trackIndex].trackRef()->chi2() / elements[trackIndex].trackRef()->ndof();
0469   const float nlost = elements[trackIndex].trackRef()->missingInnerHits();
0470   const float nLayers = elements[trackIndex].trackRef()->hitPattern().trackerLayersWithMeasurement();
0471   const float trackPt = elements[trackIndex].trackRef()->pt();
0472   const float stip = elements[trackIndex].trackRefPF()->STIP();
0473 
0474   float linkedE = 0;
0475   float linkedH = 0;
0476   std::multimap<double, unsigned int> ecalAssoTrack;
0477   block.associatedElements(
0478       trackIndex, linkData, ecalAssoTrack, reco::PFBlockElement::ECAL, reco::PFBlock::LINKTEST_ALL);
0479   std::multimap<double, unsigned int> hcalAssoTrack;
0480   block.associatedElements(
0481       trackIndex, linkData, hcalAssoTrack, reco::PFBlockElement::HCAL, reco::PFBlock::LINKTEST_ALL);
0482   if (!ecalAssoTrack.empty()) {
0483     for (auto& itecal : ecalAssoTrack) {
0484       linkedE = linkedE + elements[itecal.second].clusterRef()->energy();
0485     }
0486   }
0487   if (!hcalAssoTrack.empty()) {
0488     for (auto& ithcal : hcalAssoTrack) {
0489       linkedH = linkedH + elements[ithcal.second].clusterRef()->energy();
0490     }
0491   }
0492   const float eOverPt = linkedE / elements[trackIndex].trackRef()->pt();
0493   const float hOverPt = linkedH / elements[trackIndex].trackRef()->pt();
0494   GlobalVector rvtx(elements[trackIndex].trackRef()->innerPosition().X() - primaryVtx.x(),
0495                     elements[trackIndex].trackRef()->innerPosition().Y() - primaryVtx.y(),
0496                     elements[trackIndex].trackRef()->innerPosition().Z() - primaryVtx.z());
0497   double vtxPhi = rvtx.phi();
0498   //delta Phi between conversion vertex and track
0499   float delPhi = fabs(deltaPhi(vtxPhi, elements[trackIndex].trackRef()->innerMomentum().Phi()));
0500 
0501   float vars[] = {delPhi, nLayers, chi2, eOverPt, hOverPt, trackPt, stip, nlost};
0502 
0503   return gbrForests_.singleLeg_->GetAdaBoostClassifier(vars);
0504 }
0505 
0506 bool PFEGammaAlgo::isMuon(const reco::PFBlockElement& pfbe) {
0507   switch (pfbe.type()) {
0508     case reco::PFBlockElement::GSF: {
0509       auto& elements = _currentblock->elements();
0510       std::multimap<double, unsigned> tks;
0511       _currentblock->associatedElements(
0512           pfbe.index(), _currentlinks, tks, reco::PFBlockElement::TRACK, reco::PFBlock::LINKTEST_ALL);
0513       for (const auto& tk : tks) {
0514         if (PFMuonAlgo::isMuon(elements[tk.second])) {
0515           return true;
0516         }
0517       }
0518     } break;
0519     case reco::PFBlockElement::TRACK:
0520       return PFMuonAlgo::isMuon(pfbe);
0521       break;
0522     default:
0523       break;
0524   }
0525   return false;
0526 }
0527 
0528 PFEGammaAlgo::EgammaObjects PFEGammaAlgo::operator()(const reco::PFBlockRef& block) {
0529   LOGVERB("PFEGammaAlgo") << "Resetting PFEGammaAlgo for new block and running!" << std::endl;
0530 
0531   // candidate collections:
0532   // this starts off as an inclusive list of prototype objects built from
0533   // supercluster/ecal-driven seeds and tracker driven seeds in a block
0534   // it is then refined through by various cleanings, determining the energy
0535   // flow.
0536   // use list for constant-time removals
0537   std::list<ProtoEGObject> refinableObjects;
0538 
0539   _splayedblock.clear();
0540   _splayedblock.resize(13);  // make sure that we always have the HGCAL entry
0541 
0542   _currentblock = block;
0543   _currentlinks = block->linkData();
0544   //LOGDRESSED("PFEGammaAlgo") << *_currentblock << std::endl;
0545   LOGVERB("PFEGammaAlgo") << "Splaying block" << std::endl;
0546   //unwrap the PF block into a fast access map
0547   for (const auto& pfelement : _currentblock->elements()) {
0548     if (isMuon(pfelement))
0549       continue;  // don't allow muons in our element list
0550     if (pfelement.type() == PFBlockElement::HCAL && pfelement.clusterRef()->flags() & reco::CaloCluster::badHcalMarker)
0551       continue;  // skip also dead area markers for now
0552     const size_t itype = (size_t)pfelement.type();
0553     if (itype >= _splayedblock.size())
0554       _splayedblock.resize(itype + 1);
0555     _splayedblock[itype].emplace_back(&pfelement, true);
0556   }
0557 
0558   // show the result of splaying the tree if it's really *really* needed
0559 #ifdef PFLOW_DEBUG
0560   std::stringstream splayout;
0561   for (size_t itype = 0; itype < _splayedblock.size(); ++itype) {
0562     splayout << "\tType: " << itype << " indices: ";
0563     for (const auto& flaggedelement : _splayedblock[itype]) {
0564       splayout << flaggedelement->index() << ' ';
0565     }
0566     if (itype != _splayedblock.size() - 1)
0567       splayout << std::endl;
0568   }
0569   LOGVERB("PFEGammaAlgo") << splayout.str();
0570 #endif
0571 
0572   // precleaning of the ECAL clusters with respect to primary KF tracks
0573   // we don't allow clusters in super clusters to be locked out this way
0574   removeOrLinkECALClustersToKFTracks();
0575 
0576   initializeProtoCands(refinableObjects);
0577   LOGDRESSED("PFEGammaAlgo") << "Initialized " << refinableObjects.size() << " proto-EGamma objects" << std::endl;
0578   dumpCurrentRefinableObjects();
0579 
0580   //
0581   // now we start the refining steps
0582   //
0583   //
0584 
0585   // --- Primary Linking Step ---
0586   // since this is particle flow and we try to work from the pixels out
0587   // we start by linking the tracks together and finding the ECAL clusters
0588   for (auto& RO : refinableObjects) {
0589     // find the KF tracks associated to GSF primary tracks
0590     linkRefinableObjectGSFTracksToKFs(RO);
0591     // do the same for HCAL clusters associated to the GSF
0592     linkRefinableObjectPrimaryGSFTrackToHCAL(RO);
0593     // link secondary KF tracks associated to primary KF tracks
0594     linkRefinableObjectPrimaryKFsToSecondaryKFs(RO);
0595     // pick up clusters that are linked to the GSF primary
0596     linkRefinableObjectPrimaryGSFTrackToECAL(RO);
0597     // link associated KF to ECAL (ECAL part grabs PS clusters too if able)
0598     linkRefinableObjectKFTracksToECAL(RO);
0599     // now finally look for clusters associated to brem tangents
0600     linkRefinableObjectBremTangentsToECAL(RO);
0601   }
0602 
0603   LOGDRESSED("PFEGammaAlgo") << "Dumping after GSF and KF Track (Primary) Linking : " << std::endl;
0604   dumpCurrentRefinableObjects();
0605 
0606   // merge objects after primary linking
0607   mergeROsByAnyLink(refinableObjects);
0608 
0609   LOGDRESSED("PFEGammaAlgo") << "Dumping after first merging operation : " << std::endl;
0610   dumpCurrentRefinableObjects();
0611 
0612   // --- Secondary Linking Step ---
0613   // after this we go through the ECAL clusters on the remaining tracks
0614   // and try to link those in...
0615   for (auto& RO : refinableObjects) {
0616     // look for conversion legs
0617     linkRefinableObjectECALToSingleLegConv(RO);
0618     dumpCurrentRefinableObjects();
0619     // look for tracks that complement conversion legs
0620     linkRefinableObjectConvSecondaryKFsToSecondaryKFs(RO);
0621     // look again for ECAL clusters (this time with an e/p cut)
0622     linkRefinableObjectSecondaryKFsToECAL(RO);
0623   }
0624 
0625   LOGDRESSED("PFEGammaAlgo") << "Dumping after ECAL to Track (Secondary) Linking : " << std::endl;
0626   dumpCurrentRefinableObjects();
0627 
0628   // merge objects after primary linking
0629   mergeROsByAnyLink(refinableObjects);
0630 
0631   LOGDRESSED("PFEGammaAlgo") << "There are " << refinableObjects.size() << " after the 2nd merging step." << std::endl;
0632   dumpCurrentRefinableObjects();
0633 
0634   // -- unlinking and proto-object vetos, final sorting
0635   for (auto& RO : refinableObjects) {
0636     // remove secondary KFs (and possibly ECALs) matched to HCAL clusters
0637     unlinkRefinableObjectKFandECALMatchedToHCAL(RO, false, false);
0638     // remove secondary KFs and ECALs linked to them that have bad E/p_in
0639     // and spoil the resolution
0640     unlinkRefinableObjectKFandECALWithBadEoverP(RO);
0641     // put things back in order after partitioning
0642     std::sort(RO.ecalclusters.begin(), RO.ecalclusters.end(), [](auto const& a, auto const& b) {
0643       return (a->clusterRef()->correctedEnergy() > b->clusterRef()->correctedEnergy());
0644     });
0645     setROElectronCluster(RO);
0646   }
0647 
0648   LOGDRESSED("PFEGammaAlgo") << "There are " << refinableObjects.size() << " after the unlinking and vetos step."
0649                              << std::endl;
0650   dumpCurrentRefinableObjects();
0651 
0652   // fill the PF candidates and then build the refined SC
0653   return fillPFCandidates(refinableObjects);
0654 }
0655 
0656 void PFEGammaAlgo::initializeProtoCands(std::list<PFEGammaAlgo::ProtoEGObject>& egobjs) {
0657   // step 1: build SC based proto-candidates
0658   // in the future there will be an SC Et requirement made here to control
0659   // block size
0660   for (auto& element : _splayedblock[PFBlockElement::SC]) {
0661     LOGDRESSED("PFEGammaAlgo") << "creating SC-based proto-object" << std::endl
0662                                << "\tSC at index: " << element->index() << " has type: " << element->type()
0663                                << std::endl;
0664     element.setFlag(false);
0665     ProtoEGObject fromSC;
0666     fromSC.nBremsWithClusters = -1;
0667     fromSC.firstBrem = -1;
0668     fromSC.lateBrem = -1;
0669     fromSC.parentBlock = _currentblock;
0670     fromSC.parentSC = docast(const PFSCElement*, element.get());
0671     // splay the supercluster so we can knock out used elements
0672     bool sc_success = unwrapSuperCluster(fromSC.parentSC, fromSC.ecalclusters, fromSC.ecal2ps);
0673     if (sc_success) {
0674       /*
0675       auto ins_pos = std::lower_bound(refinableObjects.begin(),
0676                       refinableObjects.end(),
0677                       fromSC,
0678                       [&](const ProtoEGObject& a,
0679                       const ProtoEGObject& b){
0680                     const double a_en = 
0681                     a.parentSC->superClusterRef()->energy();
0682                     const double b_en = 
0683                     b.parentSC->superClusterRef()->energy();
0684                     return a_en < b_en;
0685                       });
0686       */
0687       egobjs.insert(egobjs.end(), fromSC);
0688     }
0689   }
0690   // step 2: build GSF-seed-based proto-candidates
0691   reco::GsfTrackRef gsfref_forextra;
0692   reco::TrackExtraRef gsftrk_extra;
0693   reco::ElectronSeedRef theseedref;
0694   for (auto& element : _splayedblock[PFBlockElement::GSF]) {
0695     LOGDRESSED("PFEGammaAlgo") << "creating GSF-based proto-object" << std::endl
0696                                << "\tGSF at index: " << element->index() << " has type: " << element->type()
0697                                << std::endl;
0698     const PFGSFElement* elementAsGSF = docast(const PFGSFElement*, element.get());
0699     if (elementAsGSF->trackType(reco::PFBlockElement::T_FROM_GAMMACONV)) {
0700       continue;  // for now, do not allow dedicated brems to make proto-objects
0701     }
0702     element.setFlag(false);
0703 
0704     ProtoEGObject fromGSF;
0705     fromGSF.nBremsWithClusters = -1;
0706     fromGSF.firstBrem = -1;
0707     fromGSF.lateBrem = 0;
0708     gsfref_forextra = elementAsGSF->GsftrackRef();
0709     gsftrk_extra = (gsfref_forextra.isAvailable() ? gsfref_forextra->extra() : reco::TrackExtraRef());
0710     theseedref = (gsftrk_extra.isAvailable() ? gsftrk_extra->seedRef().castTo<reco::ElectronSeedRef>()
0711                                              : reco::ElectronSeedRef());
0712     fromGSF.electronSeed = theseedref;
0713     // exception if there's no seed
0714     if (fromGSF.electronSeed.isNull() || !fromGSF.electronSeed.isAvailable()) {
0715       std::stringstream gsf_err;
0716       elementAsGSF->Dump(gsf_err, "\t");
0717       throw cms::Exception("PFEGammaAlgo::initializeProtoCands()")
0718           << "Found a GSF track with no seed! This should not happen!" << std::endl
0719           << gsf_err.str() << std::endl;
0720     }
0721     // flag this GSF element as globally used and push back the track ref
0722     // into the protocand
0723     element.setFlag(false);
0724     fromGSF.parentBlock = _currentblock;
0725     fromGSF.primaryGSFs.push_back(elementAsGSF);
0726     // add the directly matched brem tangents
0727     for (auto& brem : _splayedblock[PFBlockElement::BREM]) {
0728       float dist =
0729           _currentblock->dist(elementAsGSF->index(), brem->index(), _currentlinks, reco::PFBlock::LINKTEST_ALL);
0730       if (dist == 0.001f) {
0731         const PFBremElement* eAsBrem = docast(const PFBremElement*, brem.get());
0732         fromGSF.brems.push_back(eAsBrem);
0733         fromGSF.localMap.insert(eAsBrem, elementAsGSF);
0734         brem.setFlag(false);
0735       }
0736     }
0737     // if this track is ECAL seeded reset links or import cluster
0738     // tracker (this is pixel only, right?) driven seeds just get the GSF
0739     // track associated since this only branches for ECAL Driven seeds
0740     if (fromGSF.electronSeed->isEcalDriven()) {
0741       // step 2a: either merge with existing ProtoEG object with SC or add
0742       //          SC directly to this proto EG object if not present
0743       LOGDRESSED("PFEGammaAlgo") << "GSF-based proto-object is ECAL driven, merging SC-cand" << std::endl;
0744       LOGVERB("PFEGammaAlgo") << "ECAL Seed Ptr: " << fromGSF.electronSeed.get()
0745                               << " isAvailable: " << fromGSF.electronSeed.isAvailable()
0746                               << " isNonnull: " << fromGSF.electronSeed.isNonnull() << std::endl;
0747       SeedMatchesToProtoObject sctoseedmatch(fromGSF.electronSeed);
0748       auto objsbegin = egobjs.begin();
0749       auto objsend = egobjs.end();
0750       // this auto is a std::list<ProtoEGObject>::iterator
0751       auto clusmatch = std::find_if(objsbegin, objsend, sctoseedmatch);
0752       if (clusmatch != objsend) {
0753         fromGSF.parentSC = clusmatch->parentSC;
0754         fromGSF.ecalclusters = std::move(clusmatch->ecalclusters);
0755         fromGSF.ecal2ps = std::move(clusmatch->ecal2ps);
0756         egobjs.erase(clusmatch);
0757       } else if (fromGSF.electronSeed.isAvailable() && fromGSF.electronSeed.isNonnull()) {
0758         // link tests in the gap region can current split a gap electron
0759         // HEY THIS IS A WORK AROUND FOR A KNOWN BUG IN PFBLOCKALGO
0760         // MAYBE WE SHOULD FIX IT??????????????????????????????????
0761         LOGDRESSED("PFEGammaAlgo") << "Encountered the known GSF-SC splitting bug "
0762                                    << " in PFBlockAlgo! We should really fix this!" << std::endl;
0763       } else {  // SC was not in a earlier proto-object
0764         std::stringstream gsf_err;
0765         elementAsGSF->Dump(gsf_err, "\t");
0766         throw cms::Exception("PFEGammaAlgo::initializeProtoCands()")
0767             << "Expected SuperCluster from ECAL driven GSF seed "
0768             << "was not found in the block!" << std::endl
0769             << gsf_err.str() << std::endl;
0770       }  // supercluster in block
0771     }    // is ECAL driven seed?
0772 
0773     egobjs.insert(egobjs.end(), fromGSF);
0774   }  // end loop on GSF elements of block
0775 }
0776 
0777 bool PFEGammaAlgo::unwrapSuperCluster(const PFSCElement* thesc,
0778                                       std::vector<FlaggedPtr<const PFClusterElement>>& ecalclusters,
0779                                       ClusterMap& ecal2ps) {
0780   ecalclusters.clear();
0781   ecal2ps.clear();
0782   LOGVERB("PFEGammaAlgo") << "Pointer to SC element: 0x" << std::hex << thesc << std::dec << std::endl
0783                           << "cleared ecalclusters and ecal2ps!" << std::endl;
0784   auto ecalbegin = _splayedblock[reco::PFBlockElement::ECAL].begin();
0785   auto ecalend = _splayedblock[reco::PFBlockElement::ECAL].end();
0786   auto hgcalbegin = _splayedblock[reco::PFBlockElement::HGCAL].begin();
0787   auto hgcalend = _splayedblock[reco::PFBlockElement::HGCAL].end();
0788   if (ecalbegin == ecalend && hgcalbegin == hgcalend) {
0789     LOGERR("PFEGammaAlgo::unwrapSuperCluster()") << "There are no ECAL elements in a block with imported SC!"
0790                                                  << " This is a bug we should fix this!" << std::endl;
0791     return false;
0792   }
0793   reco::SuperClusterRef scref = thesc->superClusterRef();
0794   const bool is_pf_sc = thesc->fromPFSuperCluster();
0795   if (!(scref.isAvailable() && scref.isNonnull())) {
0796     throw cms::Exception("PFEGammaAlgo::unwrapSuperCluster()")
0797         << "SuperCluster pointed to by block element is null!" << std::endl;
0798   }
0799   LOGDRESSED("PFEGammaAlgo") << "Got a valid super cluster ref! 0x" << std::hex << scref.get() << std::dec << std::endl;
0800   const size_t nscclusters = scref->clustersSize();
0801   const size_t nscpsclusters = scref->preshowerClustersSize();
0802   size_t npfpsclusters = 0;
0803   size_t npfclusters = 0;
0804   LOGDRESSED("PFEGammaAlgo") << "Precalculated cluster multiplicities: " << nscclusters << ' ' << nscpsclusters
0805                              << std::endl;
0806   NotCloserToOther<reco::PFBlockElement::SC, reco::PFBlockElement::ECAL> ecalClustersInSC(_currentblock, thesc);
0807   NotCloserToOther<reco::PFBlockElement::SC, reco::PFBlockElement::HGCAL> hgcalClustersInSC(_currentblock, thesc);
0808   auto ecalfirstnotinsc = std::partition(ecalbegin, ecalend, ecalClustersInSC);
0809   auto hgcalfirstnotinsc = std::partition(hgcalbegin, hgcalend, hgcalClustersInSC);
0810   //reset the begin and end iterators
0811   ecalbegin = _splayedblock[reco::PFBlockElement::ECAL].begin();
0812   ecalend = _splayedblock[reco::PFBlockElement::ECAL].end();
0813 
0814   hgcalbegin = _splayedblock[reco::PFBlockElement::HGCAL].begin();
0815   hgcalend = _splayedblock[reco::PFBlockElement::HGCAL].end();
0816 
0817   //get list of associated clusters by det id and energy matching
0818   //(only needed when using non-pf supercluster)
0819   std::vector<const ClusterElement*> safePFClusters =
0820       is_pf_sc ? std::vector<const ClusterElement*>()
0821                : getSCAssociatedECALsSafe(scref, _splayedblock[reco::PFBlockElement::ECAL]);
0822 
0823   if (ecalfirstnotinsc == ecalbegin && hgcalfirstnotinsc == hgcalbegin) {
0824     LOGERR("PFEGammaAlgo::unwrapSuperCluster()") << "No associated block elements to SuperCluster!"
0825                                                  << " This is a bug we should fix!" << std::endl;
0826     return false;
0827   }
0828   npfclusters = std::distance(ecalbegin, ecalfirstnotinsc) + std::distance(hgcalbegin, hgcalfirstnotinsc);
0829   // ensure we have found the correct number of PF ecal clusters in the case
0830   // that this is a PF supercluster, otherwise all bets are off
0831   if (is_pf_sc && nscclusters != npfclusters) {
0832     std::stringstream sc_err;
0833     thesc->Dump(sc_err, "\t");
0834     throw cms::Exception("PFEGammaAlgo::unwrapSuperCluster()")
0835         << "The number of found ecal elements (" << nscclusters << ") in block is not the same as"
0836         << " the number of ecal PF clusters reported by the PFSuperCluster"
0837         << " itself (" << npfclusters << ")! This should not happen!" << std::endl
0838         << sc_err.str() << std::endl;
0839   }
0840   for (auto ecalitr = ecalbegin; ecalitr != ecalfirstnotinsc; ++ecalitr) {
0841     const PFClusterElement* elemascluster = docast(const PFClusterElement*, ecalitr->get());
0842 
0843     // reject clusters that really shouldn't be associated to the SC
0844     // (only needed when using non-pf-supercluster)
0845     if (!is_pf_sc && std::find(safePFClusters.begin(), safePFClusters.end(), elemascluster) == safePFClusters.end())
0846       continue;
0847 
0848     //add cluster
0849     ecalclusters.emplace_back(elemascluster, true);
0850     //mark cluster as used
0851     ecalitr->setFlag(false);
0852 
0853     // process the ES elements
0854     // auto is a pair<Iterator,bool> here, bool is false when placing fails
0855     auto emplaceresult = ecal2ps.emplace(elemascluster, ClusterMap::mapped_type());
0856     if (!emplaceresult.second) {
0857       std::stringstream clus_err;
0858       elemascluster->Dump(clus_err, "\t");
0859       throw cms::Exception("PFEGammaAlgo::unwrapSuperCluster()")
0860           << "List of pointers to ECAL block elements contains non-unique items!"
0861           << " This is very bad!" << std::endl
0862           << "cluster ptr = 0x" << std::hex << elemascluster << std::dec << std::endl
0863           << clus_err.str() << std::endl;
0864     }
0865     ClusterMap::mapped_type& eslist = emplaceresult.first->second;
0866     npfpsclusters += attachPSClusters(elemascluster, eslist);
0867   }  // loop over ecal elements
0868 
0869   for (auto hgcalitr = hgcalbegin; hgcalitr != hgcalfirstnotinsc; ++hgcalitr) {
0870     const PFClusterElement* elemascluster = docast(const PFClusterElement*, hgcalitr->get());
0871 
0872     // reject clusters that really shouldn't be associated to the SC
0873     // (only needed when using non-pf-supercluster)
0874     if (!is_pf_sc && std::find(safePFClusters.begin(), safePFClusters.end(), elemascluster) == safePFClusters.end())
0875       continue;
0876 
0877     //add cluster
0878     ecalclusters.emplace_back(elemascluster, true);
0879     //mark cluster as used
0880     hgcalitr->setFlag(false);
0881   }  // loop over ecal elements
0882 
0883   /*
0884    if( is_pf_sc && nscpsclusters != npfpsclusters) {
0885      std::stringstream sc_err;
0886      thesc->Dump(sc_err,"\t");
0887      throw cms::Exception("PFEGammaAlgo::unwrapSuperCluster()")
0888        << "The number of found PF preshower elements (" 
0889        << npfpsclusters << ") in block is not the same as"
0890        << " the number of preshower clusters reported by the PFSuperCluster"
0891        << " itself (" << nscpsclusters << ")! This should not happen!" 
0892        << std::endl 
0893        << sc_err.str() << std::endl;
0894    }
0895    */
0896 
0897   LOGDRESSED("PFEGammaAlgo") << " Unwrapped SC has " << npfclusters << " ECAL sub-clusters"
0898                              << " and " << npfpsclusters << " PreShower layers 1 & 2 clusters!" << std::endl;
0899   return true;
0900 }
0901 
0902 int PFEGammaAlgo::attachPSClusters(const ClusterElement* ecalclus, ClusterMap::mapped_type& eslist) {
0903   if (ecalclus->clusterRef()->layer() == PFLayer::ECAL_BARREL)
0904     return 0;
0905   edm::Ptr<reco::PFCluster> clusptr = refToPtr(ecalclus->clusterRef());
0906   EEtoPSElement ecalkey(clusptr.key(), clusptr);
0907   auto assc_ps =
0908       std::equal_range(eetops_.cbegin(), eetops_.cend(), ecalkey, [](const EEtoPSElement& a, const EEtoPSElement& b) {
0909         return a.first < b.first;
0910       });
0911   for (const auto& ps1 : _splayedblock[reco::PFBlockElement::PS1]) {
0912     edm::Ptr<reco::PFCluster> temp = refToPtr(ps1->clusterRef());
0913     for (auto pscl = assc_ps.first; pscl != assc_ps.second; ++pscl) {
0914       if (pscl->second == temp) {
0915         const ClusterElement* pstemp = docast(const ClusterElement*, ps1.get());
0916         eslist.emplace_back(pstemp);
0917       }
0918     }
0919   }
0920   for (const auto& ps2 : _splayedblock[reco::PFBlockElement::PS2]) {
0921     edm::Ptr<reco::PFCluster> temp = refToPtr(ps2->clusterRef());
0922     for (auto pscl = assc_ps.first; pscl != assc_ps.second; ++pscl) {
0923       if (pscl->second == temp) {
0924         const ClusterElement* pstemp = docast(const ClusterElement*, ps2.get());
0925         eslist.emplace_back(pstemp);
0926       }
0927     }
0928   }
0929   return eslist.size();
0930 }
0931 
0932 void PFEGammaAlgo::dumpCurrentRefinableObjects() const {
0933 #ifdef PFLOW_DEBUG
0934   edm::LogVerbatim("PFEGammaAlgo")
0935       //<< "Dumping current block: " << std::endl << *_currentblock << std::endl
0936       << "Dumping " << refinableObjects.size() << " refinable objects for this block: " << std::endl;
0937   for (const auto& ro : refinableObjects) {
0938     std::stringstream info;
0939     info << "Refinable Object:" << std::endl;
0940     if (ro.parentSC) {
0941       info << "\tSuperCluster element attached to object:" << std::endl << '\t';
0942       ro.parentSC->Dump(info, "\t");
0943       info << std::endl;
0944     }
0945     if (ro.electronSeed.isNonnull()) {
0946       info << "\tGSF element attached to object:" << std::endl;
0947       ro.primaryGSFs.front()->Dump(info, "\t");
0948       info << std::endl;
0949       info << "firstBrem : " << ro.firstBrem << " lateBrem : " << ro.lateBrem
0950            << " nBrems with cluster : " << ro.nBremsWithClusters << std::endl;
0951       ;
0952       if (ro.electronClusters.size() && ro.electronClusters[0]) {
0953         info << "electron cluster : ";
0954         ro.electronClusters[0]->Dump(info, "\t");
0955         info << std::endl;
0956       } else {
0957         info << " no electron cluster." << std::endl;
0958       }
0959     }
0960     if (ro.primaryKFs.size()) {
0961       info << "\tPrimary KF tracks attached to object: " << std::endl;
0962       for (const auto& kf : ro.primaryKFs) {
0963         kf->Dump(info, "\t");
0964         info << std::endl;
0965       }
0966     }
0967     if (ro.secondaryKFs.size()) {
0968       info << "\tSecondary KF tracks attached to object: " << std::endl;
0969       for (const auto& kf : ro.secondaryKFs) {
0970         kf->Dump(info, "\t");
0971         info << std::endl;
0972       }
0973     }
0974     if (ro.brems.size()) {
0975       info << "\tBrem tangents attached to object: " << std::endl;
0976       for (const auto& brem : ro.brems) {
0977         brem->Dump(info, "\t");
0978         info << std::endl;
0979       }
0980     }
0981     if (ro.ecalclusters.size()) {
0982       info << "\tECAL clusters attached to object: " << std::endl;
0983       for (const auto& clus : ro.ecalclusters) {
0984         clus->Dump(info, "\t");
0985         info << std::endl;
0986         if (ro.ecal2ps.find(clus) != ro.ecal2ps.end()) {
0987           for (const auto& psclus : ro.ecal2ps.at(clus)) {
0988             info << "\t\t Attached PS Cluster: ";
0989             psclus->Dump(info, "");
0990             info << std::endl;
0991           }
0992         }
0993       }
0994     }
0995     edm::LogVerbatim("PFEGammaAlgo") << info.str();
0996   }
0997 #endif
0998 }
0999 
1000 // look through our KF tracks in this block and match
1001 void PFEGammaAlgo::removeOrLinkECALClustersToKFTracks() {
1002   typedef std::multimap<double, unsigned> MatchedMap;
1003   typedef const reco::PFBlockElementGsfTrack* GsfTrackElementPtr;
1004   if (_splayedblock[reco::PFBlockElement::ECAL].empty() || _splayedblock[reco::PFBlockElement::TRACK].empty())
1005     return;
1006   MatchedMap matchedGSFs, matchedECALs;
1007   std::unordered_map<GsfTrackElementPtr, MatchedMap> gsf_ecal_cache;
1008   for (auto& kftrack : _splayedblock[reco::PFBlockElement::TRACK]) {
1009     matchedGSFs.clear();
1010     _currentblock->associatedElements(
1011         kftrack->index(), _currentlinks, matchedGSFs, reco::PFBlockElement::GSF, reco::PFBlock::LINKTEST_ALL);
1012     if (matchedGSFs.empty()) {  // only run this if we aren't associated to GSF
1013       LesserByDistance closestTrackToECAL(_currentblock, _currentlinks, kftrack.get());
1014       auto ecalbegin = _splayedblock[reco::PFBlockElement::ECAL].begin();
1015       auto ecalend = _splayedblock[reco::PFBlockElement::ECAL].end();
1016       std::partial_sort(ecalbegin, ecalbegin + 1, ecalend, closestTrackToECAL);
1017       auto& closestECAL = _splayedblock[reco::PFBlockElement::ECAL].front();
1018       const float dist =
1019           _currentblock->dist(kftrack->index(), closestECAL->index(), _currentlinks, reco::PFBlock::LINKTEST_ALL);
1020       bool inSC = false;
1021       for (auto& sc : _splayedblock[reco::PFBlockElement::SC]) {
1022         float dist_sc =
1023             _currentblock->dist(sc->index(), closestECAL->index(), _currentlinks, reco::PFBlock::LINKTEST_ALL);
1024         if (dist_sc != -1.0f) {
1025           inSC = true;
1026           break;
1027         }
1028       }
1029 
1030       if (dist != -1.0f && closestECAL.flag()) {
1031         bool gsflinked = false;
1032         // check that this cluster is not associated to a GSF track
1033         for (const auto& gsfflag : _splayedblock[reco::PFBlockElement::GSF]) {
1034           const reco::PFBlockElementGsfTrack* elemasgsf = docast(const reco::PFBlockElementGsfTrack*, gsfflag.get());
1035           if (elemasgsf->trackType(reco::PFBlockElement::T_FROM_GAMMACONV)) {
1036             continue;  // keep clusters that have a found conversion GSF near
1037           }
1038           // make sure cache exists
1039           if (!gsf_ecal_cache.count(elemasgsf)) {
1040             matchedECALs.clear();
1041             _currentblock->associatedElements(elemasgsf->index(),
1042                                               _currentlinks,
1043                                               matchedECALs,
1044                                               reco::PFBlockElement::ECAL,
1045                                               reco::PFBlock::LINKTEST_ALL);
1046             gsf_ecal_cache.emplace(elemasgsf, matchedECALs);
1047             MatchedMap().swap(matchedECALs);
1048           }
1049           const MatchedMap& ecal_matches = gsf_ecal_cache[elemasgsf];
1050           if (!ecal_matches.empty()) {
1051             if (ecal_matches.begin()->second == closestECAL->index()) {
1052               gsflinked = true;
1053               break;
1054             }
1055           }
1056         }  // loop over primary GSF tracks
1057         if (!gsflinked && !inSC) {
1058           // determine if we should remove the matched cluster
1059           const reco::PFBlockElementTrack* kfEle = docast(const reco::PFBlockElementTrack*, kftrack.get());
1060           const reco::TrackRef& trackref = kfEle->trackRef();
1061 
1062           const int nexhits = trackref->missingInnerHits();
1063           bool fromprimaryvertex = false;
1064           for (auto vtxtks = primaryVertex_.tracks_begin(); vtxtks != primaryVertex_.tracks_end(); ++vtxtks) {
1065             if (trackref == vtxtks->castTo<reco::TrackRef>()) {
1066               fromprimaryvertex = true;
1067               break;
1068             }
1069           }  // loop over tracks in primary vertex
1070              // if associated to good non-GSF matched track remove this cluster
1071           if ((PFTrackAlgoTools::isGoodForEGMPrimary(trackref->algo()) ||
1072                PFTrackAlgoTools::isGoodForEGMPrimary(trackref->originalAlgo())) &&
1073               nexhits == 0 && fromprimaryvertex) {
1074             closestECAL.setFlag(false);
1075           }
1076         }
1077       }  // found a good closest ECAL match
1078     }    // no GSF track matched to KF
1079   }      // loop over KF elements
1080 }
1081 
1082 void PFEGammaAlgo::mergeROsByAnyLink(std::list<PFEGammaAlgo::ProtoEGObject>& ROs) {
1083   if (ROs.size() < 2)
1084     return;  // nothing to do with one or zero ROs
1085   bool check_for_merge = true;
1086   while (check_for_merge) {
1087     // bugfix for early termination merging loop (15 April 2014)
1088     // check all pairwise combinations in the list
1089     // if one has a merge shuffle it to the front of the list
1090     // if there are no merges left to do we can terminate
1091     for (auto it1 = ROs.begin(); it1 != ROs.end(); ++it1) {
1092       auto find_start = it1;
1093       ++find_start;
1094       auto has_merge = std::find_if(find_start, ROs.end(), std::bind(testIfROMergableByLink, _1, *it1));
1095       if (has_merge != ROs.end() && it1 != ROs.begin()) {
1096         std::swap(*(ROs.begin()), *it1);
1097         break;
1098       }
1099     }  // ensure mergables are shuffled to the front
1100     ProtoEGObject& thefront = ROs.front();
1101     auto mergestart = ROs.begin();
1102     ++mergestart;
1103     auto nomerge = std::partition(mergestart, ROs.end(), std::bind(testIfROMergableByLink, _1, thefront));
1104     if (nomerge != mergestart) {
1105       LOGDRESSED("PFEGammaAlgo::mergeROsByAnyLink()")
1106           << "Found objects " << std::distance(mergestart, nomerge) << " to merge by links to the front!" << std::endl;
1107       for (auto roToMerge = mergestart; roToMerge != nomerge; ++roToMerge) {
1108         //bugfix! L.Gray 14 Jan 2016
1109         // -- check that the front is still mergeable!
1110         if (!thefront.ecalclusters.empty() && !roToMerge->ecalclusters.empty()) {
1111           if (thefront.ecalclusters.front()->clusterRef()->layer() !=
1112               roToMerge->ecalclusters.front()->clusterRef()->layer()) {
1113             LOGWARN("PFEGammaAlgo::mergeROsByAnyLink") << "Tried to merge EB and EE clusters! Skipping!";
1114             ROs.push_back(*roToMerge);
1115             continue;
1116           }
1117         }
1118         //end bugfix
1119         thefront.ecalclusters.insert(
1120             thefront.ecalclusters.end(), roToMerge->ecalclusters.begin(), roToMerge->ecalclusters.end());
1121         thefront.ecal2ps.insert(roToMerge->ecal2ps.begin(), roToMerge->ecal2ps.end());
1122         thefront.secondaryKFs.insert(
1123             thefront.secondaryKFs.end(), roToMerge->secondaryKFs.begin(), roToMerge->secondaryKFs.end());
1124 
1125         thefront.localMap.concatenate(roToMerge->localMap);
1126         // TO FIX -> use best (E_gsf - E_clustersum)/E_GSF
1127         if (!thefront.parentSC && roToMerge->parentSC) {
1128           thefront.parentSC = roToMerge->parentSC;
1129         }
1130         if (thefront.electronSeed.isNull() && roToMerge->electronSeed.isNonnull()) {
1131           thefront.electronSeed = roToMerge->electronSeed;
1132           thefront.primaryGSFs.insert(
1133               thefront.primaryGSFs.end(), roToMerge->primaryGSFs.begin(), roToMerge->primaryGSFs.end());
1134           thefront.primaryKFs.insert(
1135               thefront.primaryKFs.end(), roToMerge->primaryKFs.begin(), roToMerge->primaryKFs.end());
1136           thefront.brems.insert(thefront.brems.end(), roToMerge->brems.begin(), roToMerge->brems.end());
1137           thefront.electronClusters = roToMerge->electronClusters;
1138           thefront.nBremsWithClusters = roToMerge->nBremsWithClusters;
1139           thefront.firstBrem = roToMerge->firstBrem;
1140           thefront.lateBrem = roToMerge->lateBrem;
1141         } else if (thefront.electronSeed.isNonnull() && roToMerge->electronSeed.isNonnull()) {
1142           LOGDRESSED("PFEGammaAlgo::mergeROsByAnyLink")
1143               << "Need to implement proper merging of two gsf candidates!" << std::endl;
1144         }
1145       }
1146       ROs.erase(mergestart, nomerge);
1147       // put the merged element in the back of the cleaned list
1148       ROs.push_back(ROs.front());
1149       ROs.pop_front();
1150     } else {
1151       check_for_merge = false;
1152     }
1153   }
1154   LOGDRESSED("PFEGammaAlgo::mergeROsByAnyLink()")
1155       << "After merging by links there are: " << ROs.size() << " refinable EGamma objects!" << std::endl;
1156 }
1157 
1158 // pull in KF tracks associated to the RO but not closer to another
1159 // NB: in initializeProtoCands() we forced the GSF tracks not to be
1160 //     from a conversion, but we will leave a protection here just in
1161 //     case things change in the future
1162 void PFEGammaAlgo::linkRefinableObjectGSFTracksToKFs(ProtoEGObject& RO) {
1163   constexpr reco::PFBlockElement::TrackType convType = reco::PFBlockElement::T_FROM_GAMMACONV;
1164   if (_splayedblock[reco::PFBlockElement::TRACK].empty())
1165     return;
1166   auto KFbegin = _splayedblock[reco::PFBlockElement::TRACK].begin();
1167   auto KFend = _splayedblock[reco::PFBlockElement::TRACK].end();
1168   for (auto& gsfflagged : RO.primaryGSFs) {
1169     const PFGSFElement* seedtk = gsfflagged;
1170     // don't process SC-only ROs or secondary seeded ROs
1171     if (RO.electronSeed.isNull() || seedtk->trackType(convType))
1172       continue;
1173     NotCloserToOther<reco::PFBlockElement::GSF, reco::PFBlockElement::TRACK> gsfTrackToKFs(_currentblock, seedtk);
1174     // get KF tracks not closer to another and not already used
1175     auto notlinked = std::partition(KFbegin, KFend, gsfTrackToKFs);
1176     // attach tracks and set as used
1177     for (auto kft = KFbegin; kft != notlinked; ++kft) {
1178       const PFKFElement* elemaskf = docast(const PFKFElement*, kft->get());
1179       // don't care about things that aren't primaries or directly
1180       // associated secondary tracks
1181       if (isPrimaryTrack(*elemaskf, *seedtk) && !elemaskf->trackType(convType)) {
1182         kft->setFlag(false);
1183         RO.primaryKFs.push_back(elemaskf);
1184         RO.localMap.insert(seedtk, elemaskf);
1185       } else if (elemaskf->trackType(convType)) {
1186         kft->setFlag(false);
1187         RO.secondaryKFs.push_back(elemaskf);
1188         RO.localMap.insert(seedtk, elemaskf);
1189       }
1190     }  // loop on closest KFs not closer to other GSFs
1191   }    // loop on GSF primaries on RO
1192 }
1193 
1194 void PFEGammaAlgo::linkRefinableObjectPrimaryKFsToSecondaryKFs(ProtoEGObject& RO) {
1195   constexpr reco::PFBlockElement::TrackType convType = reco::PFBlockElement::T_FROM_GAMMACONV;
1196   if (_splayedblock[reco::PFBlockElement::TRACK].empty())
1197     return;
1198   auto KFbegin = _splayedblock[reco::PFBlockElement::TRACK].begin();
1199   auto KFend = _splayedblock[reco::PFBlockElement::TRACK].end();
1200   for (auto& primkf : RO.primaryKFs) {
1201     // don't process SC-only ROs or secondary seeded ROs
1202     if (primkf->trackType(convType)) {
1203       throw cms::Exception("PFEGammaAlgo::linkRefinableObjectPrimaryKFsToSecondaryKFs()")
1204           << "A KF track from conversion has been assigned as a primary!!" << std::endl;
1205     }
1206     NotCloserToOther<reco::PFBlockElement::TRACK, reco::PFBlockElement::TRACK, true> kfTrackToKFs(_currentblock,
1207                                                                                                   primkf);
1208     // get KF tracks not closer to another and not already used
1209     auto notlinked = std::partition(KFbegin, KFend, kfTrackToKFs);
1210     // attach tracks and set as used
1211     for (auto kft = KFbegin; kft != notlinked; ++kft) {
1212       const PFKFElement* elemaskf = docast(const PFKFElement*, kft->get());
1213       // don't care about things that aren't primaries or directly
1214       // associated secondary tracks
1215       if (elemaskf->trackType(convType)) {
1216         kft->setFlag(false);
1217         RO.secondaryKFs.push_back(elemaskf);
1218         RO.localMap.insert(primkf, elemaskf);
1219       }
1220     }  // loop on closest KFs not closer to other KFs
1221   }    // loop on KF primaries on RO
1222 }
1223 
1224 // try to associate the tracks to cluster elements which are not used
1225 void PFEGammaAlgo::linkRefinableObjectPrimaryGSFTrackToECAL(ProtoEGObject& RO) {
1226   if (_splayedblock[reco::PFBlockElement::ECAL].empty()) {
1227     RO.electronClusters.push_back(nullptr);
1228     return;
1229   }
1230   auto ECALbegin = _splayedblock[reco::PFBlockElement::ECAL].begin();
1231   auto ECALend = _splayedblock[reco::PFBlockElement::ECAL].end();
1232   for (auto& primgsf : RO.primaryGSFs) {
1233     NotCloserToOther<reco::PFBlockElement::GSF, reco::PFBlockElement::ECAL> gsfTracksToECALs(_currentblock, primgsf);
1234     // get set of matching ecals not already in SC
1235     auto notmatched_blk = std::partition(ECALbegin, ECALend, gsfTracksToECALs);
1236     notmatched_blk =
1237         std::partition(ECALbegin, notmatched_blk, [&primgsf](auto const& x) { return compatibleEoPOut(*x, *primgsf); });
1238     // get set of matching ecals already in the RO
1239     auto notmatched_sc = std::partition(RO.ecalclusters.begin(), RO.ecalclusters.end(), gsfTracksToECALs);
1240     notmatched_sc = std::partition(
1241         RO.ecalclusters.begin(), notmatched_sc, [&primgsf](auto const& x) { return compatibleEoPOut(*x, *primgsf); });
1242     // look inside the SC for the ECAL cluster
1243     for (auto ecal = RO.ecalclusters.begin(); ecal != notmatched_sc; ++ecal) {
1244       const PFClusterElement* elemascluster = docast(const PFClusterElement*, ecal->get());
1245       FlaggedPtr<const PFClusterElement> temp(elemascluster, true);
1246       LOGDRESSED("PFEGammaAlgo::linkGSFTracktoECAL()") << "Found a cluster already in RO by GSF extrapolation"
1247                                                        << " at ECAL surface!" << std::endl
1248                                                        << *elemascluster << std::endl;
1249 
1250       RO.localMap.insert(primgsf, temp.get());
1251     }
1252     // look outside the SC for the ecal cluster
1253     for (auto ecal = ECALbegin; ecal != notmatched_blk; ++ecal) {
1254       const PFClusterElement* elemascluster = docast(const PFClusterElement*, ecal->get());
1255       LOGDRESSED("PFEGammaAlgo::linkGSFTracktoECAL()") << "Found a cluster not already in RO by GSF extrapolation"
1256                                                        << " at ECAL surface!" << std::endl
1257                                                        << *elemascluster << std::endl;
1258       if (addPFClusterToROSafe(elemascluster, RO)) {
1259         attachPSClusters(elemascluster, RO.ecal2ps[elemascluster]);
1260         RO.localMap.insert(primgsf, elemascluster);
1261         ecal->setFlag(false);
1262       }
1263     }
1264   }
1265 }
1266 
1267 // try to associate the tracks to cluster elements which are not used
1268 void PFEGammaAlgo::linkRefinableObjectPrimaryGSFTrackToHCAL(ProtoEGObject& RO) {
1269   if (_splayedblock[reco::PFBlockElement::HCAL].empty())
1270     return;
1271   auto HCALbegin = _splayedblock[reco::PFBlockElement::HCAL].begin();
1272   auto HCALend = _splayedblock[reco::PFBlockElement::HCAL].end();
1273   for (auto& primgsf : RO.primaryGSFs) {
1274     NotCloserToOther<reco::PFBlockElement::GSF, reco::PFBlockElement::HCAL> gsfTracksToHCALs(_currentblock, primgsf);
1275     auto notmatched = std::partition(HCALbegin, HCALend, gsfTracksToHCALs);
1276     for (auto hcal = HCALbegin; hcal != notmatched; ++hcal) {
1277       const PFClusterElement* elemascluster = docast(const PFClusterElement*, hcal->get());
1278       FlaggedPtr<const PFClusterElement> temp(elemascluster, true);
1279       LOGDRESSED("PFEGammaAlgo::linkGSFTracktoECAL()")
1280           << "Found an HCAL cluster associated to GSF extrapolation" << std::endl;
1281       RO.hcalClusters.push_back(temp.get());
1282       RO.localMap.insert(primgsf, temp.get());
1283       hcal->setFlag(false);
1284     }
1285   }
1286 }
1287 
1288 // try to associate the tracks to cluster elements which are not used
1289 void PFEGammaAlgo::linkRefinableObjectKFTracksToECAL(ProtoEGObject& RO) {
1290   if (_splayedblock[reco::PFBlockElement::ECAL].empty())
1291     return;
1292   for (auto& primkf : RO.primaryKFs)
1293     linkKFTrackToECAL(primkf, RO);
1294   for (auto& secdkf : RO.secondaryKFs)
1295     linkKFTrackToECAL(secdkf, RO);
1296 }
1297 
1298 void PFEGammaAlgo::linkKFTrackToECAL(PFKFElement const* kfflagged, ProtoEGObject& RO) {
1299   std::vector<FlaggedPtr<const PFClusterElement>>& currentECAL = RO.ecalclusters;
1300   auto ECALbegin = _splayedblock[reco::PFBlockElement::ECAL].begin();
1301   auto ECALend = _splayedblock[reco::PFBlockElement::ECAL].end();
1302   NotCloserToOther<reco::PFBlockElement::TRACK, reco::PFBlockElement::ECAL> kfTrackToECALs(_currentblock, kfflagged);
1303   NotCloserToOther<reco::PFBlockElement::GSF, reco::PFBlockElement::ECAL> kfTrackGSFToECALs(_currentblock, kfflagged);
1304   //get the ECAL elements not used and not closer to another KF
1305   auto notmatched_sc = std::partition(currentECAL.begin(), currentECAL.end(), kfTrackToECALs);
1306   //get subset ECAL elements not used or closer to another GSF of any type
1307   notmatched_sc = std::partition(currentECAL.begin(), notmatched_sc, kfTrackGSFToECALs);
1308   for (auto ecalitr = currentECAL.begin(); ecalitr != notmatched_sc; ++ecalitr) {
1309     const PFClusterElement* elemascluster = docast(const PFClusterElement*, ecalitr->get());
1310     FlaggedPtr<const PFClusterElement> flaggedclus(elemascluster, true);
1311 
1312     LOGDRESSED("PFEGammaAlgo::linkKFTracktoECAL()") << "Found a cluster already in RO by KF extrapolation"
1313                                                     << " at ECAL surface!" << std::endl
1314                                                     << *elemascluster << std::endl;
1315     RO.localMap.insert(elemascluster, kfflagged);
1316   }
1317   //get the ECAL elements not used and not closer to another KF
1318   auto notmatched_blk = std::partition(ECALbegin, ECALend, kfTrackToECALs);
1319   //get subset ECAL elements not used or closer to another GSF of any type
1320   notmatched_blk = std::partition(ECALbegin, notmatched_blk, kfTrackGSFToECALs);
1321   for (auto ecalitr = ECALbegin; ecalitr != notmatched_blk; ++ecalitr) {
1322     const PFClusterElement* elemascluster = docast(const PFClusterElement*, ecalitr->get());
1323     if (addPFClusterToROSafe(elemascluster, RO)) {
1324       attachPSClusters(elemascluster, RO.ecal2ps[elemascluster]);
1325       ecalitr->setFlag(false);
1326 
1327       LOGDRESSED("PFEGammaAlgo::linkKFTracktoECAL()") << "Found a cluster not in RO by KF extrapolation"
1328                                                       << " at ECAL surface!" << std::endl
1329                                                       << *elemascluster << std::endl;
1330       RO.localMap.insert(elemascluster, kfflagged);
1331     }
1332   }
1333 }
1334 
1335 void PFEGammaAlgo::linkRefinableObjectBremTangentsToECAL(ProtoEGObject& RO) {
1336   if (RO.brems.empty())
1337     return;
1338   int FirstBrem = -1;
1339   int TrajPos = -1;
1340   int lastBremTrajPos = -1;
1341   for (auto& brem : RO.brems) {
1342     bool has_clusters = false;
1343     TrajPos = (brem->indTrajPoint()) - 2;
1344     auto ECALbegin = _splayedblock[reco::PFBlockElement::ECAL].begin();
1345     auto ECALend = _splayedblock[reco::PFBlockElement::ECAL].end();
1346     NotCloserToOther<reco::PFBlockElement::BREM, reco::PFBlockElement::ECAL> BremToECALs(_currentblock, brem);
1347     // check for late brem using clusters already in the SC
1348     auto RSCBegin = RO.ecalclusters.begin();
1349     auto RSCEnd = RO.ecalclusters.end();
1350     auto notmatched_rsc = std::partition(RSCBegin, RSCEnd, BremToECALs);
1351     for (auto ecal = RSCBegin; ecal != notmatched_rsc; ++ecal) {
1352       float deta = std::abs((*ecal)->clusterRef()->positionREP().eta() - brem->positionAtECALEntrance().eta());
1353       if (deta < 0.015) {
1354         has_clusters = true;
1355         if (lastBremTrajPos == -1 || lastBremTrajPos < TrajPos) {
1356           lastBremTrajPos = TrajPos;
1357         }
1358         if (FirstBrem == -1 || TrajPos < FirstBrem) {  // set brem information
1359           FirstBrem = TrajPos;
1360           RO.firstBrem = TrajPos;
1361         }
1362         LOGDRESSED("PFEGammaAlgo::linkBremToECAL()") << "Found a cluster already in SC linked to brem extrapolation"
1363                                                      << " at ECAL surface!" << std::endl;
1364         RO.localMap.insert(ecal->get(), brem);
1365       }
1366     }
1367     // grab new clusters from the block (ensured to not be late brem)
1368     auto notmatched_block = std::partition(ECALbegin, ECALend, BremToECALs);
1369     for (auto ecal = ECALbegin; ecal != notmatched_block; ++ecal) {
1370       float deta = std::abs((*ecal)->clusterRef()->positionREP().eta() - brem->positionAtECALEntrance().eta());
1371       if (deta < 0.015) {
1372         has_clusters = true;
1373         if (lastBremTrajPos == -1 || lastBremTrajPos < TrajPos) {
1374           lastBremTrajPos = TrajPos;
1375         }
1376         if (FirstBrem == -1 || TrajPos < FirstBrem) {  // set brem information
1377 
1378           FirstBrem = TrajPos;
1379           RO.firstBrem = TrajPos;
1380         }
1381         const PFClusterElement* elemasclus = docast(const PFClusterElement*, ecal->get());
1382         if (addPFClusterToROSafe(elemasclus, RO)) {
1383           attachPSClusters(elemasclus, RO.ecal2ps[elemasclus]);
1384 
1385           RO.localMap.insert(ecal->get(), brem);
1386           ecal->setFlag(false);
1387           LOGDRESSED("PFEGammaAlgo::linkBremToECAL()") << "Found a cluster not already associated by brem extrapolation"
1388                                                        << " at ECAL surface!" << std::endl;
1389         }
1390       }
1391     }
1392     if (has_clusters) {
1393       if (RO.nBremsWithClusters == -1)
1394         RO.nBremsWithClusters = 0;
1395       ++RO.nBremsWithClusters;
1396     }
1397   }
1398 }
1399 
1400 void PFEGammaAlgo::linkRefinableObjectConvSecondaryKFsToSecondaryKFs(ProtoEGObject& RO) {
1401   auto KFbegin = _splayedblock[reco::PFBlockElement::TRACK].begin();
1402   auto KFend = _splayedblock[reco::PFBlockElement::TRACK].end();
1403   auto BeginROskfs = RO.secondaryKFs.begin();
1404   auto EndROskfs = RO.secondaryKFs.end();
1405   auto ronotconv = std::partition(BeginROskfs, EndROskfs, [](auto const& x) { return x->trackType(ConvType); });
1406   size_t convkfs_end = std::distance(BeginROskfs, ronotconv);
1407   for (size_t idx = 0; idx < convkfs_end; ++idx) {
1408     auto const& secKFs =
1409         RO.secondaryKFs;  //we want the entry at the index but we allocate to secondaryKFs in loop which invalidates all iterators, references and pointers, hence we need to get the entry fresh each time
1410     NotCloserToOther<reco::PFBlockElement::TRACK, reco::PFBlockElement::TRACK, true> TracksToTracks(_currentblock,
1411                                                                                                     secKFs[idx]);
1412     auto notmatched = std::partition(KFbegin, KFend, TracksToTracks);
1413     notmatched = std::partition(KFbegin, notmatched, [](auto const& x) { return x->trackType(ConvType); });
1414     for (auto kf = KFbegin; kf != notmatched; ++kf) {
1415       const reco::PFBlockElementTrack* elemaskf = docast(const reco::PFBlockElementTrack*, kf->get());
1416       RO.secondaryKFs.push_back(elemaskf);
1417       RO.localMap.insert(secKFs[idx], kf->get());
1418       kf->setFlag(false);
1419     }
1420   }
1421 }
1422 
1423 void PFEGammaAlgo::linkRefinableObjectECALToSingleLegConv(ProtoEGObject& RO) {
1424   auto KFbegin = _splayedblock[reco::PFBlockElement::TRACK].begin();
1425   auto KFend = _splayedblock[reco::PFBlockElement::TRACK].end();
1426   for (auto& ecal : RO.ecalclusters) {
1427     NotCloserToOther<reco::PFBlockElement::ECAL, reco::PFBlockElement::TRACK, true> ECALToTracks(_currentblock,
1428                                                                                                  ecal.get());
1429     auto notmatchedkf = std::partition(KFbegin, KFend, ECALToTracks);
1430     auto notconvkf = std::partition(KFbegin, notmatchedkf, [](auto const& x) { return x->trackType(ConvType); });
1431     // add identified KF conversion tracks
1432     for (auto kf = KFbegin; kf != notconvkf; ++kf) {
1433       const reco::PFBlockElementTrack* elemaskf = docast(const reco::PFBlockElementTrack*, kf->get());
1434       RO.secondaryKFs.push_back(elemaskf);
1435       RO.localMap.insert(ecal.get(), elemaskf);
1436       kf->setFlag(false);
1437     }
1438     // go through non-conv-identified kfs and check MVA to add conversions
1439     for (auto kf = notconvkf; kf != notmatchedkf; ++kf) {
1440       float mvaval = evaluateSingleLegMVA(_currentblock, primaryVertex_, (*kf)->index());
1441       if (mvaval > cfg_.mvaConvCut) {
1442         const reco::PFBlockElementTrack* elemaskf = docast(const reco::PFBlockElementTrack*, kf->get());
1443         RO.secondaryKFs.push_back(elemaskf);
1444         RO.localMap.insert(ecal.get(), elemaskf);
1445         kf->setFlag(false);
1446 
1447         RO.singleLegConversionMvaMap.emplace(elemaskf, mvaval);
1448       }
1449     }
1450   }
1451 }
1452 
1453 void PFEGammaAlgo::linkRefinableObjectSecondaryKFsToECAL(ProtoEGObject& RO) {
1454   auto ECALbegin = _splayedblock[reco::PFBlockElement::ECAL].begin();
1455   auto ECALend = _splayedblock[reco::PFBlockElement::ECAL].end();
1456   for (auto& skf : RO.secondaryKFs) {
1457     NotCloserToOther<reco::PFBlockElement::TRACK, reco::PFBlockElement::ECAL, false> TracksToECALwithCut(
1458         _currentblock, skf, 1.5f);
1459     auto notmatched = std::partition(ECALbegin, ECALend, TracksToECALwithCut);
1460     for (auto ecal = ECALbegin; ecal != notmatched; ++ecal) {
1461       const reco::PFBlockElementCluster* elemascluster = docast(const reco::PFBlockElementCluster*, ecal->get());
1462       if (addPFClusterToROSafe(elemascluster, RO)) {
1463         attachPSClusters(elemascluster, RO.ecal2ps[elemascluster]);
1464         RO.localMap.insert(skf, elemascluster);
1465         ecal->setFlag(false);
1466       }
1467     }
1468   }
1469 }
1470 
1471 PFEGammaAlgo::EgammaObjects PFEGammaAlgo::fillPFCandidates(const std::list<PFEGammaAlgo::ProtoEGObject>& ROs) {
1472   EgammaObjects output;
1473 
1474   // reserve output collections
1475   output.candidates.reserve(ROs.size());
1476   output.candidateExtras.reserve(ROs.size());
1477   output.refinedSuperClusters.reserve(ROs.size());
1478 
1479   for (auto& RO : ROs) {
1480     if (RO.ecalclusters.empty() && !cfg_.produceEGCandsWithNoSuperCluster)
1481       continue;
1482 
1483     reco::PFCandidate cand;
1484     reco::PFCandidateEGammaExtra xtra;
1485     if (!RO.primaryGSFs.empty() || !RO.primaryKFs.empty()) {
1486       cand.setPdgId(-11);  // anything with a primary track is an electron
1487     } else {
1488       cand.setPdgId(22);  // anything with no primary track is a photon
1489     }
1490     if (!RO.primaryKFs.empty()) {
1491       cand.setCharge(RO.primaryKFs[0]->trackRef()->charge());
1492       xtra.setKfTrackRef(RO.primaryKFs[0]->trackRef());
1493       cand.setTrackRef(RO.primaryKFs[0]->trackRef());
1494       cand.setVertex(RO.primaryKFs[0]->trackRef()->vertex());
1495       cand.addElementInBlock(_currentblock, RO.primaryKFs[0]->index());
1496     }
1497     if (!RO.primaryGSFs.empty()) {
1498       cand.setCharge(RO.primaryGSFs[0]->GsftrackRef()->chargeMode());
1499       xtra.setGsfTrackRef(RO.primaryGSFs[0]->GsftrackRef());
1500       cand.setGsfTrackRef(RO.primaryGSFs[0]->GsftrackRef());
1501       cand.setVertex(RO.primaryGSFs[0]->GsftrackRef()->vertex());
1502       cand.addElementInBlock(_currentblock, RO.primaryGSFs[0]->index());
1503     }
1504     if (RO.parentSC) {
1505       xtra.setSuperClusterPFECALRef(RO.parentSC->superClusterRef());
1506       // we'll set to the refined supercluster back up in the producer
1507       cand.setSuperClusterRef(RO.parentSC->superClusterRef());
1508       xtra.setSuperClusterRef(RO.parentSC->superClusterRef());
1509       cand.addElementInBlock(_currentblock, RO.parentSC->index());
1510     }
1511     // add brems
1512     for (const auto& brem : RO.brems) {
1513       cand.addElementInBlock(_currentblock, brem->index());
1514     }
1515     // add clusters and ps clusters
1516     for (const auto& ecal : RO.ecalclusters) {
1517       const PFClusterElement* clus = ecal.get();
1518       cand.addElementInBlock(_currentblock, clus->index());
1519       if (RO.ecal2ps.count(clus)) {
1520         for (auto& psclus : RO.ecal2ps.at(clus)) {
1521           cand.addElementInBlock(_currentblock, psclus->index());
1522         }
1523       }
1524     }
1525     // add secondary tracks
1526     for (const auto& kf : RO.secondaryKFs) {
1527       cand.addElementInBlock(_currentblock, kf->index());
1528       const reco::ConversionRefVector& convrefs = kf->convRefs();
1529       bool no_conv_ref = true;
1530       for (const auto& convref : convrefs) {
1531         if (convref.isNonnull() && convref.isAvailable()) {
1532           xtra.addConversionRef(convref);
1533           no_conv_ref = false;
1534         }
1535       }
1536       if (no_conv_ref) {
1537         //single leg conversions
1538 
1539         //look for stored mva value in map or else recompute
1540         const auto& mvavalmapped = RO.singleLegConversionMvaMap.find(kf);
1541         //FIXME: Abuse single mva value to store both provenance and single leg mva score
1542         //by storing 3.0 + mvaval
1543         float mvaval = (mvavalmapped != RO.singleLegConversionMvaMap.end()
1544                             ? mvavalmapped->second
1545                             : 3.0 + evaluateSingleLegMVA(_currentblock, primaryVertex_, kf->index()));
1546 
1547         xtra.addSingleLegConvTrackRefMva(std::make_pair(kf->trackRef(), mvaval));
1548       }
1549     }
1550 
1551     // build the refined supercluster from those clusters left in the cand
1552     output.refinedSuperClusters.push_back(buildRefinedSuperCluster(RO));
1553 
1554     //*TODO* cluster time is not reliable at the moment, so only use track timing
1555     float trkTime = 0, trkTimeErr = -1;
1556     if (!RO.primaryGSFs.empty() && RO.primaryGSFs[0]->isTimeValid()) {
1557       trkTime = RO.primaryGSFs[0]->time();
1558       trkTimeErr = RO.primaryGSFs[0]->timeError();
1559     } else if (!RO.primaryKFs.empty() && RO.primaryKFs[0]->isTimeValid()) {
1560       trkTime = RO.primaryKFs[0]->time();
1561       trkTimeErr = RO.primaryKFs[0]->timeError();
1562     }
1563     if (trkTimeErr >= 0) {
1564       cand.setTime(trkTime, trkTimeErr);
1565     }
1566 
1567     const reco::SuperCluster& the_sc = output.refinedSuperClusters.back();
1568     // with the refined SC in hand we build a naive candidate p4
1569     // and set the candidate ECAL position to either the barycenter of the
1570     // supercluster (if super-cluster present) or the seed of the
1571     // new SC generated by the EGAlgo
1572     const double scE = the_sc.energy();
1573     if (scE != 0.0) {
1574       const math::XYZPoint& seedPos = the_sc.seed()->position();
1575       math::XYZVector egDir = the_sc.position() - primaryVertex_.position();
1576       egDir = egDir.Unit();
1577       cand.setP4(math::XYZTLorentzVector(scE * egDir.x(), scE * egDir.y(), scE * egDir.z(), scE));
1578       math::XYZPointF ecalPOS_f(seedPos.x(), seedPos.y(), seedPos.z());
1579       cand.setPositionAtECALEntrance(ecalPOS_f);
1580       cand.setEcalEnergy(the_sc.rawEnergy(), the_sc.energy());
1581     } else if (cfg_.produceEGCandsWithNoSuperCluster && !RO.primaryGSFs.empty()) {
1582       const PFGSFElement* gsf = RO.primaryGSFs[0];
1583       const reco::GsfTrackRef& gref = gsf->GsftrackRef();
1584       math::XYZTLorentzVector p4(gref->pxMode(), gref->pyMode(), gref->pzMode(), gref->pMode());
1585       cand.setP4(p4);
1586       cand.setPositionAtECALEntrance(gsf->positionAtECALEntrance());
1587     } else if (cfg_.produceEGCandsWithNoSuperCluster && !RO.primaryKFs.empty()) {
1588       const PFKFElement* kf = RO.primaryKFs[0];
1589       reco::TrackRef kref = RO.primaryKFs[0]->trackRef();
1590       math::XYZTLorentzVector p4(kref->px(), kref->py(), kref->pz(), kref->p());
1591       cand.setP4(p4);
1592       cand.setPositionAtECALEntrance(kf->positionAtECALEntrance());
1593     }
1594     const float eleMVAValue = calculateEleMVA(RO, xtra);
1595     fillExtraInfo(RO, xtra);
1596     //std::cout << "PFEG eleMVA: " << eleMVAValue << std::endl;
1597     xtra.setMVA(eleMVAValue);
1598     cand.set_mva_e_pi(eleMVAValue);
1599     output.candidates.push_back(cand);
1600     output.candidateExtras.push_back(xtra);
1601   }
1602 
1603   return output;
1604 }
1605 
1606 float PFEGammaAlgo::calculateEleMVA(const PFEGammaAlgo::ProtoEGObject& ro, reco::PFCandidateEGammaExtra& xtra) const {
1607   if (ro.primaryGSFs.empty()) {
1608     return -2.0f;
1609   }
1610   const PFGSFElement* gsfElement = ro.primaryGSFs.front();
1611   const PFKFElement* kfElement = nullptr;
1612   if (!ro.primaryKFs.empty()) {
1613     kfElement = ro.primaryKFs.front();
1614   }
1615   auto const& refGsf = gsfElement->GsftrackRef();
1616   reco::TrackRef refKf;
1617   constexpr float mEl = 0.000511;
1618   const double eInGsf = std::hypot(refGsf->pMode(), mEl);
1619   double dEtGsfEcal = 1e6;
1620   double sigmaEtaEta = 1e-14;
1621   const double eneHcalGsf =
1622       std::accumulate(ro.hcalClusters.begin(), ro.hcalClusters.end(), 0.0, [](const double a, auto const& b) {
1623         return a + b->clusterRef()->energy();
1624       });
1625   if (!ro.primaryKFs.empty()) {
1626     refKf = ro.primaryKFs.front()->trackRef();
1627   }
1628   const double eOutGsf = gsfElement->Pout().t();
1629   const double etaOutGsf = gsfElement->positionAtECALEntrance().eta();
1630   double firstEcalGsfEnergy{0.0};
1631   double otherEcalGsfEnergy{0.0};
1632   double ecalBremEnergy{0.0};
1633   //shower shape of cluster closest to gsf track
1634   std::vector<const reco::PFCluster*> gsfCluster;
1635   for (const auto& ecal : ro.ecalclusters) {
1636     const double cenergy = ecal->clusterRef()->correctedEnergy();
1637     bool hasgsf = ro.localMap.contains(gsfElement, ecal.get());
1638     bool haskf = ro.localMap.contains(kfElement, ecal.get());
1639     bool hasbrem = false;
1640     for (const auto& brem : ro.brems) {
1641       if (ro.localMap.contains(brem, ecal.get())) {
1642         hasbrem = true;
1643       }
1644     }
1645     if (hasbrem && ecal.get() != ro.electronClusters[0]) {
1646       ecalBremEnergy += cenergy;
1647     }
1648     if (!hasbrem && ecal.get() != ro.electronClusters[0]) {
1649       if (hasgsf)
1650         otherEcalGsfEnergy += cenergy;
1651       if (haskf)
1652         ecalBremEnergy += cenergy;  // from conv. brem!
1653       if (!(hasgsf || haskf))
1654         otherEcalGsfEnergy += cenergy;  // stuff from SC
1655     }
1656   }
1657 
1658   if (ro.electronClusters[0]) {
1659     reco::PFClusterRef cref = ro.electronClusters[0]->clusterRef();
1660     xtra.setGsfElectronClusterRef(_currentblock, *(ro.electronClusters[0]));
1661     firstEcalGsfEnergy = cref->correctedEnergy();
1662     dEtGsfEcal = cref->positionREP().eta() - etaOutGsf;
1663     gsfCluster.push_back(cref.get());
1664     PFClusterWidthAlgo pfwidth(gsfCluster);
1665     sigmaEtaEta = pfwidth.pflowSigmaEtaEta();
1666   }
1667 
1668   // brem sequence information
1669   float firstBrem{-1.0f};
1670   float earlyBrem{-1.0f};
1671   float lateBrem{-1.0f};
1672   if (ro.nBremsWithClusters > 0) {
1673     firstBrem = ro.firstBrem;
1674     earlyBrem = ro.firstBrem < 4 ? 1.0f : 0.0f;
1675     lateBrem = ro.lateBrem == 1 ? 1.0f : 0.0f;
1676   }
1677   xtra.setEarlyBrem(earlyBrem);
1678   xtra.setLateBrem(lateBrem);
1679   if (firstEcalGsfEnergy > 0.0) {
1680     if (refGsf.isNonnull()) {
1681       xtra.setGsfTrackPout(gsfElement->Pout());
1682       // normalization observables
1683       const float ptGsf = refGsf->ptMode();
1684       const float etaGsf = refGsf->etaMode();
1685       // tracking observables
1686       const double ptModeErrorGsf = refGsf->ptModeError();
1687       float ptModeErrOverPtGsf = (ptModeErrorGsf > 0. ? ptModeErrorGsf / ptGsf : 1.0);
1688       float chi2Gsf = refGsf->normalizedChi2();
1689       float dPtOverPtGsf = (ptGsf - gsfElement->Pout().pt()) / ptGsf;
1690       // kalman filter vars
1691       float nHitKf = refKf.isNonnull() ? refKf->hitPattern().trackerLayersWithMeasurement() : 0;
1692       float chi2Kf = refKf.isNonnull() ? refKf->normalizedChi2() : -0.01;
1693 
1694       //tracker + calorimetry observables
1695       float eTotPinMode = (firstEcalGsfEnergy + otherEcalGsfEnergy + ecalBremEnergy) / eInGsf;
1696       float eGsfPoutMode = firstEcalGsfEnergy / eOutGsf;
1697       float eTotBremPinPoutMode = (ecalBremEnergy + otherEcalGsfEnergy) / (eInGsf - eOutGsf);
1698       float dEtaGsfEcalClust = std::abs(dEtGsfEcal);
1699       float logSigmaEtaEta = std::log(sigmaEtaEta);
1700       float hOverHe = eneHcalGsf / (eneHcalGsf + firstEcalGsfEnergy);
1701 
1702       xtra.setDeltaEta(dEtaGsfEcalClust);
1703       xtra.setSigmaEtaEta(sigmaEtaEta);
1704       xtra.setHadEnergy(eneHcalGsf);
1705 
1706       // Apply bounds to variables and calculate MVA
1707       dPtOverPtGsf = std::clamp(dPtOverPtGsf, -0.2f, 1.0f);
1708       ptModeErrOverPtGsf = std::min(ptModeErrOverPtGsf, 0.3f);
1709       chi2Gsf = std::min(chi2Gsf, 10.0f);
1710       chi2Kf = std::min(chi2Kf, 10.0f);
1711       eTotPinMode = std::clamp(eTotPinMode, 0.0f, 5.0f);
1712       eGsfPoutMode = std::clamp(eGsfPoutMode, 0.0f, 5.0f);
1713       eTotBremPinPoutMode = std::clamp(eTotBremPinPoutMode, 0.0f, 5.0f);
1714       dEtaGsfEcalClust = std::min(dEtaGsfEcalClust, 0.1f);
1715       logSigmaEtaEta = std::max(logSigmaEtaEta, -14.0f);
1716 
1717       // not used for moment, weird behavior of variable
1718       //float dPtOverPtKf = refKf.isNonnull() ? (refKf->pt() - refKf->outerPt())/refKf->pt() : -0.01;
1719       //dPtOverPtKf       = std::clamp(dPtOverPtKf,-0.2f, 1.0f);
1720 
1721       /*
1722  *      To be used for debugging:
1723  *      pretty-print the PFEgamma electron MVA input variables
1724  *
1725  *      std::cout << " **** PFEG BDT observables ****" << endl;
1726  *      std::cout << " < Normalization > " << endl;
1727  *      std::cout << " ptGsf " << ptGsf << " Pin " << eInGsf
1728  *        << " Pout " << eOutGsf << " etaGsf " << etaGsf << endl;
1729  *      std::cout << " < PureTracking > " << endl;
1730  *      std::cout << " ptModeErrOverPtGsf " << ptModeErrOverPtGsf 
1731  *        << " dPtOverPtGsf " << dPtOverPtGsf
1732  *        << " chi2Gsf " << chi2Gsf
1733  *        << " nhit_gsf " << nhit_gsf
1734  *        << " dPtOverPtKf " << dPtOverPtKf
1735  *        << " chi2Kf " << chi2Kf 
1736  *        << " nHitKf " << nHitKf <<  endl;
1737  *      std::cout << " < track-ecal-hcal-ps " << endl;
1738  *      std::cout << " eTotPinMode " << eTotPinMode 
1739  *        << " eGsfPoutMode " << eGsfPoutMode
1740  *        << " eTotBremPinPoutMode " << eTotBremPinPoutMode
1741  *        << " dEtaGsfEcalClust " << dEtaGsfEcalClust 
1742  *        << " logSigmaEtaEta " << logSigmaEtaEta
1743  *        << " hOverHe " << hOverHe << " Hcal energy " << eneHcalGsf
1744  *        << " lateBrem " << lateBrem
1745  *        << " firstBrem " << firstBrem << endl;
1746  */
1747 
1748       float vars[] = {std::log(ptGsf),
1749                       etaGsf,
1750                       ptModeErrOverPtGsf,
1751                       dPtOverPtGsf,
1752                       chi2Gsf,
1753                       nHitKf,
1754                       chi2Kf,
1755                       eTotPinMode,
1756                       eGsfPoutMode,
1757                       eTotBremPinPoutMode,
1758                       dEtaGsfEcalClust,
1759                       logSigmaEtaEta,
1760                       hOverHe,
1761                       lateBrem,
1762                       firstBrem};
1763 
1764       return gbrForests_.ele_->GetAdaBoostClassifier(vars);
1765     }
1766   }
1767   return -2.0f;
1768 }
1769 
1770 void PFEGammaAlgo::fillExtraInfo(const ProtoEGObject& RO, reco::PFCandidateEGammaExtra& xtra) {
1771   // add tracks associated to clusters that are not T_FROM_GAMMACONV
1772   // info about single-leg convs is already save, so just veto in loops
1773   auto KFbegin = _splayedblock[reco::PFBlockElement::TRACK].begin();
1774   auto KFend = _splayedblock[reco::PFBlockElement::TRACK].end();
1775   for (auto& ecal : RO.ecalclusters) {
1776     NotCloserToOther<reco::PFBlockElement::ECAL, reco::PFBlockElement::TRACK, true> ECALToTracks(_currentblock,
1777                                                                                                  ecal.get());
1778     auto notmatchedkf = std::partition(KFbegin, KFend, ECALToTracks);
1779     auto notconvkf = std::partition(KFbegin, notmatchedkf, [](auto const& x) { return x->trackType(ConvType); });
1780     // go through non-conv-identified kfs and check MVA to add conversions
1781     for (auto kf = notconvkf; kf != notmatchedkf; ++kf) {
1782       const reco::PFBlockElementTrack* elemaskf = docast(const reco::PFBlockElementTrack*, kf->get());
1783       xtra.addExtraNonConvTrack(_currentblock, *elemaskf);
1784     }
1785   }
1786 }
1787 
1788 // currently stolen from PFECALSuperClusterAlgo, we should
1789 // try to factor this correctly since the operation is the same in
1790 // both places...
1791 reco::SuperCluster PFEGammaAlgo::buildRefinedSuperCluster(const PFEGammaAlgo::ProtoEGObject& RO) {
1792   if (RO.ecalclusters.empty()) {
1793     return reco::SuperCluster(0.0, math::XYZPoint(0, 0, 0));
1794   }
1795 
1796   bool isEE = false;
1797   // need the vector of raw pointers for a PF width class
1798   std::vector<const reco::PFCluster*> bare_ptrs;
1799   // calculate necessary parameters and build the SC
1800   double posX(0), posY(0), posZ(0), rawSCEnergy(0), corrSCEnergy(0), corrPSEnergy(0), ps1_energy(0.0), ps2_energy(0.0);
1801   for (auto& clus : RO.ecalclusters) {
1802     double ePS1 = 0;
1803     double ePS2 = 0;
1804     isEE = PFLayer::ECAL_ENDCAP == clus->clusterRef()->layer();
1805     auto clusptr = edm::refToPtr<reco::PFClusterCollection>(clus->clusterRef());
1806     bare_ptrs.push_back(clusptr.get());
1807 
1808     const double cluseraw = clusptr->energy();
1809     double cluscalibe = clusptr->correctedEnergy();
1810     const math::XYZPoint& cluspos = clusptr->position();
1811     posX += cluseraw * cluspos.X();
1812     posY += cluseraw * cluspos.Y();
1813     posZ += cluseraw * cluspos.Z();
1814     // update EE calibrated super cluster energies
1815     if (isEE && RO.ecal2ps.count(clus.get())) {
1816       const auto& psclusters = RO.ecal2ps.at(clus.get());
1817 
1818       std::vector<reco::PFCluster const*> psClusterPointers;
1819       psClusterPointers.reserve(psclusters.size());
1820       for (auto const& psc : psclusters) {
1821         psClusterPointers.push_back(psc->clusterRef().get());
1822       }
1823       auto calibratedEnergies = thePFEnergyCalibration_.calibrateEndcapClusterEnergies(
1824           *clusptr, psClusterPointers, channelStatus_, cfg_.applyCrackCorrections);
1825       cluscalibe = calibratedEnergies.clusterEnergy;
1826       ePS1 = calibratedEnergies.ps1Energy;
1827       ePS2 = calibratedEnergies.ps2Energy;
1828     }
1829     if (ePS1 == -1.)
1830       ePS1 = 0;
1831     if (ePS2 == -1.)
1832       ePS2 = 0;
1833 
1834     rawSCEnergy += cluseraw;
1835     corrSCEnergy += cluscalibe;
1836     ps1_energy += ePS1;
1837     ps2_energy += ePS2;
1838     corrPSEnergy += ePS1 + ePS2;
1839   }
1840   posX /= rawSCEnergy;
1841   posY /= rawSCEnergy;
1842   posZ /= rawSCEnergy;
1843 
1844   // now build the supercluster
1845   reco::SuperCluster new_sc(corrSCEnergy, math::XYZPoint(posX, posY, posZ));
1846 
1847   auto clusptr = edm::refToPtr<reco::PFClusterCollection>(RO.ecalclusters.front()->clusterRef());
1848   new_sc.setCorrectedEnergy(corrSCEnergy);
1849   new_sc.setSeed(clusptr);
1850   new_sc.setPreshowerEnergyPlane1(ps1_energy);
1851   new_sc.setPreshowerEnergyPlane2(ps2_energy);
1852   new_sc.setPreshowerEnergy(corrPSEnergy);
1853   for (const auto& clus : RO.ecalclusters) {
1854     clusptr = edm::refToPtr<reco::PFClusterCollection>(clus->clusterRef());
1855     new_sc.addCluster(clusptr);
1856     auto& hits_and_fractions = clusptr->hitsAndFractions();
1857     for (auto& hit_and_fraction : hits_and_fractions) {
1858       new_sc.addHitAndFraction(hit_and_fraction.first, hit_and_fraction.second);
1859     }
1860     // put the preshower stuff back in later
1861     if (RO.ecal2ps.count(clus.get())) {
1862       const auto& cluspsassociation = RO.ecal2ps.at(clus.get());
1863       // EE rechits should be uniquely matched to sets of pre-shower
1864       // clusters at this point, so we throw an exception if otherwise
1865       // now wrapped in EDM debug flags
1866       for (const auto& pscluselem : cluspsassociation) {
1867         edm::Ptr<reco::PFCluster> psclus = edm::refToPtr<reco::PFClusterCollection>(pscluselem->clusterRef());
1868 #ifdef PFFLOW_DEBUG
1869         auto found_pscluster =
1870             std::find(new_sc.preshowerClustersBegin(), new_sc.preshowerClustersEnd(), reco::CaloClusterPtr(psclus));
1871         if (found_pscluster == new_sc.preshowerClustersEnd()) {
1872 #endif
1873           new_sc.addPreshowerCluster(psclus);
1874 #ifdef PFFLOW_DEBUG
1875         } else {
1876           throw cms::Exception("PFECALSuperClusterAlgo::buildSuperCluster")
1877               << "Found a PS cluster matched to more than one EE cluster!" << std::endl
1878               << std::hex << psclus.get() << " == " << found_pscluster->get() << std::dec << std::endl;
1879         }
1880 #endif
1881       }
1882     }
1883   }
1884 
1885   // calculate linearly weighted cluster widths
1886   PFClusterWidthAlgo pfwidth(bare_ptrs);
1887   new_sc.setEtaWidth(pfwidth.pflowEtaWidth());
1888   new_sc.setPhiWidth(pfwidth.pflowPhiWidth());
1889 
1890   // cache the value of the raw energy
1891   new_sc.rawEnergy();
1892 
1893   return new_sc;
1894 }
1895 
1896 void PFEGammaAlgo::unlinkRefinableObjectKFandECALWithBadEoverP(ProtoEGObject& RO) {
1897   // this only means something for ROs with a primary GSF track
1898   if (RO.primaryGSFs.empty())
1899     return;
1900   // need energy sums to tell if we've added crap or not
1901   const double Pin_gsf = RO.primaryGSFs.front()->GsftrackRef()->pMode();
1902   const double gsfOuterEta = RO.primaryGSFs.front()->positionAtECALEntrance().Eta();
1903   double tot_ecal = 0.0;
1904   std::vector<double> min_brem_dists;
1905   std::vector<double> closest_brem_eta;
1906   // first get the total ecal energy (we should replace this with a cache)
1907   for (const auto& ecal : RO.ecalclusters) {
1908     tot_ecal += ecal->clusterRef()->correctedEnergy();
1909     // we also need to look at the minimum distance to brems
1910     // since energetic brems will be closer to the brem than the track
1911     double min_brem_dist = 5000.0;
1912     double eta = -999.0;
1913     for (const auto& brem : RO.brems) {
1914       const float dist = _currentblock->dist(brem->index(), ecal->index(), _currentlinks, reco::PFBlock::LINKTEST_ALL);
1915       if (dist < min_brem_dist && dist != -1.0f) {
1916         min_brem_dist = dist;
1917         eta = brem->positionAtECALEntrance().Eta();
1918       }
1919     }
1920     min_brem_dists.push_back(min_brem_dist);
1921     closest_brem_eta.push_back(eta);
1922   }
1923 
1924   // loop through the ECAL clusters and remove ECAL clusters matched to
1925   // secondary track either in *or* out of the SC if the E/pin is bad
1926   for (auto secd_kf = RO.secondaryKFs.begin(); secd_kf != RO.secondaryKFs.end(); ++secd_kf) {
1927     reco::TrackRef trkRef = (*secd_kf)->trackRef();
1928     const float secpin = (*secd_kf)->trackRef()->p();
1929     bool remove_this_kf = false;
1930     for (auto ecal = RO.ecalclusters.begin(); ecal != RO.ecalclusters.end(); ++ecal) {
1931       size_t bremidx = std::distance(RO.ecalclusters.begin(), ecal);
1932       const float minbremdist = min_brem_dists[bremidx];
1933       const double ecalenergy = (*ecal)->clusterRef()->correctedEnergy();
1934       const double Epin = ecalenergy / secpin;
1935       const double detaGsf = std::abs(gsfOuterEta - (*ecal)->clusterRef()->positionREP().Eta());
1936       const double detaBrem = std::abs(closest_brem_eta[bremidx] - (*ecal)->clusterRef()->positionREP().Eta());
1937 
1938       bool kf_matched = RO.localMap.contains(ecal->get(), *secd_kf);
1939 
1940       const float tkdist =
1941           _currentblock->dist((*secd_kf)->index(), (*ecal)->index(), _currentlinks, reco::PFBlock::LINKTEST_ALL);
1942 
1943       // do not reject this track if it is closer to a brem than the
1944       // secondary track, or if it lies in the delta-eta plane with the
1945       // gsf track or if it is in the dEta plane with the brems
1946       if (Epin > 3 && kf_matched && tkdist != -1.0f && tkdist < minbremdist && detaGsf > 0.05 && detaBrem > 0.015) {
1947         double res_with = std::abs((tot_ecal - Pin_gsf) / Pin_gsf);
1948         double res_without = std::abs((tot_ecal - ecalenergy - Pin_gsf) / Pin_gsf);
1949         if (res_without < res_with) {
1950           LOGDRESSED("PFEGammaAlgo") << " REJECTED_RES totenergy " << tot_ecal << " Pin_gsf " << Pin_gsf
1951                                      << " cluster to secondary " << ecalenergy << " res_with " << res_with
1952                                      << " res_without " << res_without << std::endl;
1953           tot_ecal -= ecalenergy;
1954           remove_this_kf = true;
1955           ecal = RO.ecalclusters.erase(ecal);
1956           if (ecal == RO.ecalclusters.end())
1957             break;
1958         }
1959       }
1960     }
1961     if (remove_this_kf) {
1962       secd_kf = RO.secondaryKFs.erase(secd_kf);
1963       if (secd_kf == RO.secondaryKFs.end())
1964         break;
1965     }
1966   }
1967 }
1968 
1969 void PFEGammaAlgo::unlinkRefinableObjectKFandECALMatchedToHCAL(ProtoEGObject& RO,
1970                                                                bool removeFreeECAL,
1971                                                                bool removeSCEcal) {
1972   std::vector<bool> cluster_in_sc;
1973   auto ecal_begin = RO.ecalclusters.begin();
1974   auto ecal_end = RO.ecalclusters.end();
1975   auto hcal_begin = _splayedblock[reco::PFBlockElement::HCAL].begin();
1976   auto hcal_end = _splayedblock[reco::PFBlockElement::HCAL].end();
1977   for (auto secd_kf = RO.secondaryKFs.begin(); secd_kf != RO.secondaryKFs.end(); ++secd_kf) {
1978     bool remove_this_kf = false;
1979     NotCloserToOther<reco::PFBlockElement::TRACK, reco::PFBlockElement::HCAL> tracksToHCALs(_currentblock, *secd_kf);
1980     reco::TrackRef trkRef = (*secd_kf)->trackRef();
1981 
1982     bool goodTrack = PFTrackAlgoTools::isGoodForEGM(trkRef->algo());
1983     const float secpin = trkRef->p();
1984 
1985     for (auto ecal = ecal_begin; ecal != ecal_end; ++ecal) {
1986       const double ecalenergy = (*ecal)->clusterRef()->correctedEnergy();
1987       // first check if the cluster is in the SC (use dist calc for fastness)
1988       const size_t clus_idx = std::distance(ecal_begin, ecal);
1989       if (cluster_in_sc.size() < clus_idx + 1) {
1990         float dist = -1.0f;
1991         if (RO.parentSC) {
1992           dist = _currentblock->dist((*secd_kf)->index(), (*ecal)->index(), _currentlinks, reco::PFBlock::LINKTEST_ALL);
1993         }
1994         cluster_in_sc.push_back(dist != -1.0f);
1995       }
1996 
1997       // if we've found a secondary KF that matches this ecal cluster
1998       // now we see if it is matched to HCAL
1999       // if it is matched to an HCAL cluster we take different
2000       // actions if the cluster was in an SC or not
2001       if (RO.localMap.contains(ecal->get(), *secd_kf)) {
2002         auto hcal_matched = std::partition(hcal_begin, hcal_end, tracksToHCALs);
2003         for (auto hcalclus = hcal_begin; hcalclus != hcal_matched; ++hcalclus) {
2004           const reco::PFBlockElementCluster* clusthcal = docast(const reco::PFBlockElementCluster*, hcalclus->get());
2005           const double hcalenergy = clusthcal->clusterRef()->energy();
2006           const double hpluse = ecalenergy + hcalenergy;
2007           const bool isHoHE = ((hcalenergy / hpluse) > 0.1 && goodTrack);
2008           const bool isHoE = (hcalenergy > ecalenergy);
2009           const bool isPoHE = (secpin > hpluse);
2010           if (cluster_in_sc[clus_idx]) {
2011             if (isHoE || isPoHE) {
2012               LOGDRESSED("PFEGammaAlgo") << "REJECTED TRACK FOR H/E or P/(H+E), CLUSTER IN SC"
2013                                          << " H/H+E " << (hcalenergy / hpluse) << " H/E " << (hcalenergy > ecalenergy)
2014                                          << " P/(H+E) " << (secpin / hpluse) << " HCAL ENE " << hcalenergy
2015                                          << " ECAL ENE " << ecalenergy << " secPIN " << secpin << " Algo Track "
2016                                          << trkRef->algo() << std::endl;
2017               remove_this_kf = true;
2018             }
2019           } else {
2020             if (isHoHE) {
2021               LOGDRESSED("PFEGammaAlgo") << "REJECTED TRACK FOR H/H+E, CLUSTER NOT IN SC"
2022                                          << " H/H+E " << (hcalenergy / hpluse) << " H/E " << (hcalenergy > ecalenergy)
2023                                          << " P/(H+E) " << (secpin / hpluse) << " HCAL ENE " << hcalenergy
2024                                          << " ECAL ENE " << ecalenergy << " secPIN " << secpin << " Algo Track "
2025                                          << trkRef->algo() << std::endl;
2026               remove_this_kf = true;
2027             }
2028           }
2029         }
2030       }
2031     }
2032     if (remove_this_kf) {
2033       secd_kf = RO.secondaryKFs.erase(secd_kf);
2034       if (secd_kf == RO.secondaryKFs.end())
2035         break;
2036     }
2037   }
2038 }
2039 
2040 bool PFEGammaAlgo::isPrimaryTrack(const reco::PFBlockElementTrack& KfEl, const reco::PFBlockElementGsfTrack& GsfEl) {
2041   bool isPrimary = false;
2042 
2043   const GsfPFRecTrackRef& gsfPfRef = GsfEl.GsftrackRefPF();
2044 
2045   if (gsfPfRef.isNonnull()) {
2046     const PFRecTrackRef& kfPfRef = KfEl.trackRefPF();
2047     PFRecTrackRef kfPfRef_fromGsf = (*gsfPfRef).kfPFRecTrackRef();
2048     if (kfPfRef.isNonnull() && kfPfRef_fromGsf.isNonnull()) {
2049       reco::TrackRef kfref = (*kfPfRef).trackRef();
2050       reco::TrackRef kfref_fromGsf = (*kfPfRef_fromGsf).trackRef();
2051       if (kfref.isNonnull() && kfref_fromGsf.isNonnull()) {
2052         if (kfref == kfref_fromGsf)
2053           isPrimary = true;
2054       }
2055     }
2056   }
2057 
2058   return isPrimary;
2059 }