CMS 3D CMS Logo
PFEGammaAlgo.cc
Go to the documentation of this file.
1 #include "RecoParticleFlow/PFProducer/interface/PFEGammaAlgo.h"
2 #include "RecoParticleFlow/PFProducer/interface/PFMuonAlgo.h"
3 #include "DataFormats/ParticleFlowReco/interface/PFRecHitFraction.h"
4 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
5 #include "RecoParticleFlow/PFClusterTools/interface/ClusterClusterMapping.h"
6 #include "DataFormats/ParticleFlowReco/interface/PFLayer.h"
7 #include "DataFormats/ParticleFlowReco/interface/PFRecHit.h"
8 #include "DataFormats/EcalDetId/interface/EBDetId.h"
9 #include "DataFormats/EcalDetId/interface/EEDetId.h"
10 #include "DataFormats/EcalDetId/interface/ESDetId.h"
11 #include "DataFormats/EgammaCandidates/interface/Photon.h"
12 #include "RecoParticleFlow/PFClusterTools/interface/PFEnergyCalibration.h"
13 #include "RecoParticleFlow/PFClusterTools/interface/PFEnergyResolution.h"
14 #include "RecoParticleFlow/PFClusterTools/interface/PFClusterWidthAlgo.h"
15 #include "RecoParticleFlow/PFTracking/interface/PFTrackAlgoTools.h"
16 #include "DataFormats/Common/interface/RefToPtr.h"
17 #include "RecoEcal/EgammaCoreTools/interface/Mustache.h"
18 #include "DataFormats/Math/interface/deltaPhi.h"
19 #include "DataFormats/Math/interface/deltaR.h"
20 
21 #include <TFile.h>
22 #include <TVector2.h>
23 #include <iomanip>
24 #include <algorithm>
25 #include <functional>
26 #include <numeric>
27 #include <TMath.h>
28 
29 // include combinations header (was never incorporated in boost)
30 #include "combination.hpp"
31 
32 // just for now do this
33 //#define PFLOW_DEBUG
34 
35 #ifdef PFLOW_DEBUG
36 #define docast(x, y) dynamic_cast<x>(y)
37 #define LOGVERB(x) edm::LogVerbatim(x)
38 #define LOGWARN(x) edm::LogWarning(x)
39 #define LOGERR(x) edm::LogError(x)
40 #define LOGDRESSED(x) edm::LogInfo(x)
41 #else
42 #define docast(x, y) reinterpret_cast<x>(y)
43 #define LOGVERB(x) LogTrace(x)
44 #define LOGWARN(x) edm::LogWarning(x)
45 #define LOGERR(x) edm::LogError(x)
46 #define LOGDRESSED(x) LogDebug(x)
47 #endif
48 
49 using namespace std;
50 using namespace reco;
51 using namespace std::placeholders; // for _1, _2, _3...
52 
53 namespace {
54  typedef PFEGammaAlgo::PFSCElement SCElement;
55  typedef PFEGammaAlgo::EEtoPSAssociation EEtoPSAssociation;
56  typedef std::pair<CaloClusterPtr::key_type, CaloClusterPtr> EEtoPSElement;
57  typedef PFEGammaAlgo::PFClusterElement ClusterElement;
58 
59  class SeedMatchesToProtoObject {
60  public:
61  SeedMatchesToProtoObject(const reco::ElectronSeedRef& s)
62  : scfromseed_(s->caloCluster().castTo<reco::SuperClusterRef>()) {
63  ispfsc_ = false;
64  if (scfromseed_.isNonnull()) {
65  const edm::Ptr<reco::PFCluster> testCast(scfromseed_->seed());
66  ispfsc_ = testCast.isNonnull();
67  }
68  }
69  bool operator()(const PFEGammaAlgo::ProtoEGObject& po) {
70  if (scfromseed_.isNull() || !po.parentSC)
71  return false;
72  if (ispfsc_) {
73  return (scfromseed_->seed() == po.parentSC->superClusterRef()->seed());
74  }
75  return (scfromseed_->seed()->seed() == po.parentSC->superClusterRef()->seed()->seed());
76  }
77 
78  private:
79  const reco::SuperClusterRef scfromseed_;
80  bool ispfsc_;
81  };
82 
83  template <bool useConvs = false>
84  bool elementNotCloserToOther(const reco::PFBlockRef& block,
85  const PFBlockElement::Type& keytype,
86  const size_t key,
87  const PFBlockElement::Type& valtype,
88  const size_t test,
89  const float EoPin_cut = 1.0e6) {
90  constexpr reco::PFBlockElement::TrackType ConvType = reco::PFBlockElement::T_FROM_GAMMACONV;
91  // this is inside out but I just want something that works right now
92  switch (keytype) {
93  case reco::PFBlockElement::GSF: {
94  const reco::PFBlockElementGsfTrack* elemasgsf =
95  docast(const reco::PFBlockElementGsfTrack*, &(block->elements()[key]));
96  if (elemasgsf && valtype == PFBlockElement::ECAL) {
97  const ClusterElement* elemasclus = reinterpret_cast<const ClusterElement*>(&(block->elements()[test]));
98  float cluster_e = elemasclus->clusterRef()->correctedEnergy();
99  float trk_pin = elemasgsf->Pin().P();
100  if (cluster_e / trk_pin > EoPin_cut) {
101  LOGDRESSED("elementNotCloserToOther") << "GSF track failed EoP cut to match with cluster!";
102  return false;
103  }
104  }
105  } break;
106  case reco::PFBlockElement::TRACK: {
107  const reco::PFBlockElementTrack* elemaskf = docast(const reco::PFBlockElementTrack*, &(block->elements()[key]));
108  if (elemaskf && valtype == PFBlockElement::ECAL) {
109  const ClusterElement* elemasclus = reinterpret_cast<const ClusterElement*>(&(block->elements()[test]));
110  float cluster_e = elemasclus->clusterRef()->correctedEnergy();
111  float trk_pin = std::sqrt(elemaskf->trackRef()->innerMomentum().mag2());
112  if (cluster_e / trk_pin > EoPin_cut) {
113  LOGDRESSED("elementNotCloserToOther") << "KF track failed EoP cut to match with cluster!";
114  return false;
115  }
116  }
117  } break;
118  default:
119  break;
120  }
121 
122  const float dist = block->dist(key, test, block->linkData(), reco::PFBlock::LINKTEST_ALL);
123  if (dist == -1.0f)
124  return false; // don't associate non-linked elems
125  std::multimap<double, unsigned> dists_to_val;
126  block->associatedElements(test, block->linkData(), dists_to_val, keytype, reco::PFBlock::LINKTEST_ALL);
127 
128  for (const auto& valdist : dists_to_val) {
129  const size_t idx = valdist.second;
130  // check track types for conversion info
131  switch (keytype) {
132  case reco::PFBlockElement::GSF: {
133  const reco::PFBlockElementGsfTrack* elemasgsf =
134  docast(const reco::PFBlockElementGsfTrack*, &(block->elements()[idx]));
135  if (!useConvs && elemasgsf->trackType(ConvType))
136  return false;
137  if (elemasgsf && valtype == PFBlockElement::ECAL) {
138  const ClusterElement* elemasclus = docast(const ClusterElement*, &(block->elements()[test]));
139  float cluster_e = elemasclus->clusterRef()->correctedEnergy();
140  float trk_pin = elemasgsf->Pin().P();
141  if (cluster_e / trk_pin > EoPin_cut)
142  continue;
143  }
144  } break;
145  case reco::PFBlockElement::TRACK: {
146  const reco::PFBlockElementTrack* elemaskf =
147  docast(const reco::PFBlockElementTrack*, &(block->elements()[idx]));
148  if (!useConvs && elemaskf->trackType(ConvType))
149  return false;
150  if (elemaskf && valtype == PFBlockElement::ECAL) {
151  const ClusterElement* elemasclus = reinterpret_cast<const ClusterElement*>(&(block->elements()[test]));
152  float cluster_e = elemasclus->clusterRef()->correctedEnergy();
153  float trk_pin = std::sqrt(elemaskf->trackRef()->innerMomentum().mag2());
154  if (cluster_e / trk_pin > EoPin_cut)
155  continue;
156  }
157  } break;
158  default:
159  break;
160  }
161  if (valdist.first < dist && idx != key) {
162  LOGDRESSED("elementNotCloserToOther")
163  << "key element of type " << keytype << " is closer to another element of type" << valtype << std::endl;
164  return false; // false if closer element of specified type found
165  }
166  }
167  return true;
168  }
169 
170  template <class Element1, class Element2>
171  bool compatibleEoPOut(const Element1& e, const Element2& comp) {
172  if (PFBlockElement::ECAL != e.type()) {
173  return false;
174  }
175  const ClusterElement& elemascluster = docast(ClusterElement const&, e);
176  const float gsf_eta_diff = std::abs(comp.positionAtECALEntrance().eta() - comp.Pout().eta());
177  const reco::PFClusterRef& cRef = elemascluster.clusterRef();
178  return (gsf_eta_diff <= 0.3 && cRef->energy() / comp.Pout().t() <= 5);
179  }
180 
181  constexpr reco::PFBlockElement::TrackType ConvType = reco::PFBlockElement::T_FROM_GAMMACONV;
182 
183  template <PFBlockElement::Type keytype, PFBlockElement::Type valtype, bool useConv = false>
184 
185  struct NotCloserToOther {
186  const reco::PFBlockElement* comp;
187  const reco::PFBlockRef& block;
188  const float EoPin_cut;
189  NotCloserToOther(const reco::PFBlockRef& b, const reco::PFBlockElement* e, const float EoPcut = 1.0e6)
190  : comp(e), block(b), EoPin_cut(EoPcut) {}
191  template <class T>
192  bool operator()(const T& e) {
193  if (!e.flag() || valtype != e->type())
194  return false;
195  return elementNotCloserToOther<useConv>(block, keytype, comp->index(), valtype, e->index(), EoPin_cut);
196  }
197  };
198 
199  struct LesserByDistance {
200  const reco::PFBlockElement* comp;
201  const reco::PFBlockRef& block;
202  const reco::PFBlock::LinkData& links;
203  LesserByDistance(const reco::PFBlockRef& b, const reco::PFBlock::LinkData& l, const reco::PFBlockElement* e)
204  : comp(e), block(b), links(l) {}
205  bool operator()(FlaggedPtr<const reco::PFBlockElement> const& e1,
206  FlaggedPtr<const reco::PFBlockElement> const& e2) {
207  double dist1 = block->dist(comp->index(), e1->index(), links, reco::PFBlock::LINKTEST_ALL);
208  double dist2 = block->dist(comp->index(), e2->index(), links, reco::PFBlock::LINKTEST_ALL);
209  dist1 = (dist1 == -1.0 ? 1e6 : dist1);
210  dist2 = (dist2 == -1.0 ? 1e6 : dist2);
211  return dist1 < dist2;
212  }
213  };
214 
215  bool isROLinkedByClusterOrTrack(const PFEGammaAlgo::ProtoEGObject& RO1, const PFEGammaAlgo::ProtoEGObject& RO2) {
216  // also don't allow ROs where both have clusters
217  // and GSF tracks to merge (10 Dec 2013)
218  if (!RO1.primaryGSFs.empty() && !RO2.primaryGSFs.empty()) {
219  LOGDRESSED("isROLinkedByClusterOrTrack") << "cannot merge, both have GSFs!" << std::endl;
220  return false;
221  }
222  // don't allow EB/EE to mix (11 Sept 2013)
223  if (!RO1.ecalclusters.empty() && !RO2.ecalclusters.empty()) {
224  if (RO1.ecalclusters.front()->clusterRef()->layer() != RO2.ecalclusters.front()->clusterRef()->layer()) {
225  LOGDRESSED("isROLinkedByClusterOrTrack") << "cannot merge, different ECAL types!" << std::endl;
226  return false;
227  }
228  }
229  const reco::PFBlockRef& blk = RO1.parentBlock;
230  bool not_closer;
231  // check links track -> cluster
232  for (const auto& cluster : RO1.ecalclusters) {
233  for (const auto& primgsf : RO2.primaryGSFs) {
234  not_closer = elementNotCloserToOther(blk, cluster->type(), cluster->index(), primgsf->type(), primgsf->index());
235  if (not_closer) {
236  LOGDRESSED("isROLinkedByClusterOrTrack") << "merged by cluster to primary GSF" << std::endl;
237  return true;
238  } else {
239  LOGDRESSED("isROLinkedByClusterOrTrack") << "cluster to primary GSF failed since"
240  << " cluster closer to another GSF" << std::endl;
241  }
242  }
243  for (const auto& primkf : RO2.primaryKFs) {
244  not_closer = elementNotCloserToOther(blk, cluster->type(), cluster->index(), primkf->type(), primkf->index());
245  if (not_closer) {
246  LOGDRESSED("isROLinkedByClusterOrTrack") << "merged by cluster to primary KF" << std::endl;
247  return true;
248  }
249  }
250  for (const auto& secdkf : RO2.secondaryKFs) {
251  not_closer = elementNotCloserToOther(blk, cluster->type(), cluster->index(), secdkf->type(), secdkf->index());
252  if (not_closer) {
253  LOGDRESSED("isROLinkedByClusterOrTrack") << "merged by cluster to secondary KF" << std::endl;
254  return true;
255  }
256  }
257  // check links brem -> cluster
258  for (const auto& brem : RO2.brems) {
259  not_closer = elementNotCloserToOther(blk, cluster->type(), cluster->index(), brem->type(), brem->index());
260  if (not_closer) {
261  LOGDRESSED("isROLinkedByClusterOrTrack") << "merged by cluster to brem KF" << std::endl;
262  return true;
263  }
264  }
265  }
266  // check links primary gsf -> secondary kf
267  for (const auto& primgsf : RO1.primaryGSFs) {
268  for (const auto& secdkf : RO2.secondaryKFs) {
269  not_closer = elementNotCloserToOther(blk, primgsf->type(), primgsf->index(), secdkf->type(), secdkf->index());
270  if (not_closer) {
271  LOGDRESSED("isROLinkedByClusterOrTrack") << "merged by GSF to secondary KF" << std::endl;
272  return true;
273  }
274  }
275  }
276  // check links primary kf -> secondary kf
277  for (const auto& primkf : RO1.primaryKFs) {
278  for (const auto& secdkf : RO2.secondaryKFs) {
279  not_closer = elementNotCloserToOther(blk, primkf->type(), primkf->index(), secdkf->type(), secdkf->index());
280  if (not_closer) {
281  LOGDRESSED("isROLinkedByClusterOrTrack") << "merged by primary KF to secondary KF" << std::endl;
282  return true;
283  }
284  }
285  }
286  // check links secondary kf -> secondary kf
287  for (const auto& secdkf1 : RO1.secondaryKFs) {
288  for (const auto& secdkf2 : RO2.secondaryKFs) {
289  not_closer =
290  elementNotCloserToOther<true>(blk, secdkf1->type(), secdkf1->index(), secdkf2->type(), secdkf2->index());
291  if (not_closer) {
292  LOGDRESSED("isROLinkedByClusterOrTrack") << "merged by secondary KF to secondary KF" << std::endl;
293  return true;
294  }
295  }
296  }
297  return false;
298  }
299 
300  bool testIfROMergableByLink(const PFEGammaAlgo::ProtoEGObject& ro, PFEGammaAlgo::ProtoEGObject& comp) {
301  const bool result = (isROLinkedByClusterOrTrack(comp, ro) || isROLinkedByClusterOrTrack(ro, comp));
302  return result;
303  }
304 
305  std::vector<const ClusterElement*> getSCAssociatedECALsSafe(
306  const reco::SuperClusterRef& scref, std::vector<FlaggedPtr<const reco::PFBlockElement>>& ecals) {
307  std::vector<const ClusterElement*> cluster_list;
308  auto sccl = scref->clustersBegin();
309  auto scend = scref->clustersEnd();
310  auto pfc = ecals.begin();
311  auto pfcend = ecals.end();
312  for (; sccl != scend; ++sccl) {
313  std::vector<const ClusterElement*> matched_pfcs;
314  const double eg_energy = (*sccl)->energy();
315 
316  for (pfc = ecals.begin(); pfc != pfcend; ++pfc) {
317  const ClusterElement* pfcel = docast(const ClusterElement*, pfc->get());
318  const bool matched = ClusterClusterMapping::overlap(**sccl, *(pfcel->clusterRef()));
319  // need to protect against high energy clusters being attached
320  // to low-energy SCs
321  if (matched && pfcel->clusterRef()->energy() < 1.2 * scref->energy()) {
322  matched_pfcs.push_back(pfcel);
323  }
324  }
325  std::sort(matched_pfcs.begin(), matched_pfcs.end());
326 
327  double min_residual = 1e6;
328  std::vector<const ClusterElement*> best_comb;
329  for (size_t i = 1; i <= matched_pfcs.size(); ++i) {
330  //now we find the combination of PF-clusters which
331  //has the smallest energy residual with respect to the
332  //EG-cluster we are looking at now
333  do {
334  double energy = std::accumulate(
335  matched_pfcs.begin(), matched_pfcs.begin() + i, 0.0, [](const double a, const ClusterElement* c) {
336  return a + c->clusterRef()->energy();
337  });
338  const double resid = std::abs(energy - eg_energy);
339  if (resid < min_residual) {
340  best_comb.clear();
341  best_comb.reserve(i);
342  min_residual = resid;
343  best_comb.insert(best_comb.begin(), matched_pfcs.begin(), matched_pfcs.begin() + i);
344  }
345  } while (notboost::next_combination(matched_pfcs.begin(), matched_pfcs.begin() + i, matched_pfcs.end()));
346  }
347  for (const auto& clelem : best_comb) {
348  if (std::find(cluster_list.begin(), cluster_list.end(), clelem) == cluster_list.end()) {
349  cluster_list.push_back(clelem);
350  }
351  }
352  }
353  return cluster_list;
354  }
355  bool addPFClusterToROSafe(const ClusterElement* cl, PFEGammaAlgo::ProtoEGObject& RO) {
356  if (RO.ecalclusters.empty()) {
357  RO.ecalclusters.emplace_back(cl, true);
358  return true;
359  } else {
360  const PFLayer::Layer clayer = cl->clusterRef()->layer();
361  const PFLayer::Layer blayer = RO.ecalclusters.back()->clusterRef()->layer();
362  if (clayer == blayer) {
363  RO.ecalclusters.emplace_back(cl, true);
364  return true;
365  }
366  }
367  return false;
368  }
369 
370  // sets the cluster best associated to the GSF track
371  // leave it null if no GSF track
372  void setROElectronCluster(PFEGammaAlgo::ProtoEGObject& RO) {
373  if (RO.ecalclusters.empty())
374  return;
375  RO.lateBrem = -1;
376  RO.firstBrem = -1;
377  RO.nBremsWithClusters = -1;
378  const reco::PFBlockElementBrem *firstBrem = nullptr, *lastBrem = nullptr;
379  const reco::PFBlockElementCluster *bremCluster = nullptr, *gsfCluster = nullptr, *kfCluster = nullptr,
380  *gsfCluster_noassc = nullptr;
381  const reco::PFBlockRef& parent = RO.parentBlock;
382  int nBremClusters = 0;
383  constexpr float maxDist = 1e6;
384  float mDist_gsf(maxDist), mDist_gsf_noassc(maxDist), mDist_kf(maxDist);
385  for (const auto& cluster : RO.ecalclusters) {
386  for (const auto& gsf : RO.primaryGSFs) {
387  const bool hasclu =
388  elementNotCloserToOther(parent, gsf->type(), gsf->index(), cluster->type(), cluster->index());
389  const float deta = std::abs(cluster->clusterRef()->positionREP().eta() - gsf->positionAtECALEntrance().eta());
390  const float dphi = std::abs(
391  TVector2::Phi_mpi_pi(cluster->clusterRef()->positionREP().phi() - gsf->positionAtECALEntrance().phi()));
392  const float dist = std::hypot(deta, dphi);
393  if (hasclu && dist < mDist_gsf) {
394  gsfCluster = cluster.get();
395  mDist_gsf = dist;
396  } else if (dist < mDist_gsf_noassc) {
397  gsfCluster_noassc = cluster.get();
398  mDist_gsf_noassc = dist;
399  }
400  }
401  for (const auto& kf : RO.primaryKFs) {
402  const bool hasclu = elementNotCloserToOther(parent, kf->type(), kf->index(), cluster->type(), cluster->index());
403  const float dist = parent->dist(cluster->index(), kf->index(), parent->linkData(), reco::PFBlock::LINKTEST_ALL);
404  if (hasclu && dist < mDist_kf) {
405  kfCluster = cluster.get();
406  mDist_kf = dist;
407  }
408  }
409  for (const auto& brem : RO.brems) {
410  const bool hasclu =
411  elementNotCloserToOther(parent, brem->type(), brem->index(), cluster->type(), cluster->index());
412  if (hasclu) {
413  ++nBremClusters;
414  if (!firstBrem || (firstBrem->indTrajPoint() - 2 > brem->indTrajPoint() - 2)) {
415  firstBrem = brem;
416  }
417  if (!lastBrem || (lastBrem->indTrajPoint() - 2 < brem->indTrajPoint() - 2)) {
418  lastBrem = brem;
419  bremCluster = cluster.get();
420  }
421  }
422  }
423  }
424  if (!gsfCluster && !kfCluster && !bremCluster) {
425  gsfCluster = gsfCluster_noassc;
426  }
427  RO.nBremsWithClusters = nBremClusters;
428  RO.lateBrem = 0;
429  if (gsfCluster) {
430  RO.electronClusters.push_back(gsfCluster);
431  } else if (kfCluster) {
432  RO.electronClusters.push_back(kfCluster);
433  }
434  if (bremCluster && !gsfCluster && !kfCluster) {
435  RO.electronClusters.push_back(bremCluster);
436  }
437  if (firstBrem && RO.ecalclusters.size() > 1) {
438  RO.firstBrem = firstBrem->indTrajPoint() - 2;
439  if (bremCluster == gsfCluster)
440  RO.lateBrem = 1;
441  }
442  }
443 } // namespace
444 
445 PFEGammaAlgo::PFEGammaAlgo(const PFEGammaAlgo::PFEGConfigInfo& cfg,
446  GBRForests const& gbrForests,
447  EEtoPSAssociation const& eetops,
448  ESEEIntercalibConstants const& esEEInterCalib,
449  ESChannelStatus const& channelStatus,
450  reco::Vertex const& primaryVertex)
451  : gbrForests_(gbrForests),
452  eetops_(eetops),
453  cfg_(cfg),
454  primaryVertex_(primaryVertex),
455  channelStatus_(channelStatus) {
456  thePFEnergyCalibration_.initAlphaGamma_ESplanes_fromDB(&esEEInterCalib);
457 }
458 
459 float PFEGammaAlgo::evaluateSingleLegMVA(const reco::PFBlockRef& blockRef,
460  const reco::Vertex& primaryVtx,
461  unsigned int trackIndex) {
462  const reco::PFBlock& block = *blockRef;
463  const edm::OwnVector<reco::PFBlockElement>& elements = block.elements();
464  //use this to store linkdata in the associatedElements function below
465  const PFBlock::LinkData& linkData = block.linkData();
466  //calculate MVA Variables
467  const float chi2 = elements[trackIndex].trackRef()->chi2() / elements[trackIndex].trackRef()->ndof();
468  const float nlost = elements[trackIndex].trackRef()->hitPattern().numberOfLostHits(HitPattern::MISSING_INNER_HITS);
469  const float nLayers = elements[trackIndex].trackRef()->hitPattern().trackerLayersWithMeasurement();
470  const float trackPt = elements[trackIndex].trackRef()->pt();
471  const float stip = elements[trackIndex].trackRefPF()->STIP();
472 
473  float linkedE = 0;
474  float linkedH = 0;
475  std::multimap<double, unsigned int> ecalAssoTrack;
476  block.associatedElements(
477  trackIndex, linkData, ecalAssoTrack, reco::PFBlockElement::ECAL, reco::PFBlock::LINKTEST_ALL);
478  std::multimap<double, unsigned int> hcalAssoTrack;
479  block.associatedElements(
480  trackIndex, linkData, hcalAssoTrack, reco::PFBlockElement::HCAL, reco::PFBlock::LINKTEST_ALL);
481  if (!ecalAssoTrack.empty()) {
482  for (auto& itecal : ecalAssoTrack) {
483  linkedE = linkedE + elements[itecal.second].clusterRef()->energy();
484  }
485  }
486  if (!hcalAssoTrack.empty()) {
487  for (auto& ithcal : hcalAssoTrack) {
488  linkedH = linkedH + elements[ithcal.second].clusterRef()->energy();
489  }
490  }
491  const float eOverPt = linkedE / elements[trackIndex].trackRef()->pt();
492  const float hOverPt = linkedH / elements[trackIndex].trackRef()->pt();
493  GlobalVector rvtx(elements[trackIndex].trackRef()->innerPosition().X() - primaryVtx.x(),
494  elements[trackIndex].trackRef()->innerPosition().Y() - primaryVtx.y(),
495  elements[trackIndex].trackRef()->innerPosition().Z() - primaryVtx.z());
496  double vtxPhi = rvtx.phi();
497  //delta Phi between conversion vertex and track
498  float delPhi = fabs(deltaPhi(vtxPhi, elements[trackIndex].trackRef()->innerMomentum().Phi()));
499 
500  float vars[] = {delPhi, nLayers, chi2, eOverPt, hOverPt, trackPt, stip, nlost};
501 
502  return gbrForests_.singleLeg_->GetAdaBoostClassifier(vars);
503 }
504 
505 bool PFEGammaAlgo::isMuon(const reco::PFBlockElement& pfbe) {
506  switch (pfbe.type()) {
507  case reco::PFBlockElement::GSF: {
508  auto& elements = _currentblock->elements();
509  std::multimap<double, unsigned> tks;
510  _currentblock->associatedElements(
511  pfbe.index(), _currentlinks, tks, reco::PFBlockElement::TRACK, reco::PFBlock::LINKTEST_ALL);
512  for (const auto& tk : tks) {
513  if (PFMuonAlgo::isMuon(elements[tk.second])) {
514  return true;
515  }
516  }
517  } break;
518  case reco::PFBlockElement::TRACK:
519  return PFMuonAlgo::isMuon(pfbe);
520  break;
521  default:
522  break;
523  }
524  return false;
525 }
526 
527 PFEGammaAlgo::EgammaObjects PFEGammaAlgo::operator()(const reco::PFBlockRef& block) {
528  LOGVERB("PFEGammaAlgo") << "Resetting PFEGammaAlgo for new block and running!" << std::endl;
529 
530  // candidate collections:
531  // this starts off as an inclusive list of prototype objects built from
532  // supercluster/ecal-driven seeds and tracker driven seeds in a block
533  // it is then refined through by various cleanings, determining the energy
534  // flow.
535  // use list for constant-time removals
536  std::list<ProtoEGObject> refinableObjects;
537 
538  _splayedblock.clear();
539  _splayedblock.resize(13); // make sure that we always have the HGCAL entry
540 
541  _currentblock = block;
542  _currentlinks = block->linkData();
543  //LOGDRESSED("PFEGammaAlgo") << *_currentblock << std::endl;
544  LOGVERB("PFEGammaAlgo") << "Splaying block" << std::endl;
545  //unwrap the PF block into a fast access map
546  for (const auto& pfelement : _currentblock->elements()) {
547  if (isMuon(pfelement))
548  continue; // don't allow muons in our element list
549  if (pfelement.type() == PFBlockElement::HCAL && pfelement.clusterRef()->flags() & reco::CaloCluster::badHcalMarker)
550  continue; // skip also dead area markers for now
551  const size_t itype = (size_t)pfelement.type();
552  if (itype >= _splayedblock.size())
553  _splayedblock.resize(itype + 1);
554  _splayedblock[itype].emplace_back(&pfelement, true);
555  }
556 
557  // show the result of splaying the tree if it's really *really* needed
558 #ifdef PFLOW_DEBUG
559  std::stringstream splayout;
560  for (size_t itype = 0; itype < _splayedblock.size(); ++itype) {
561  splayout << "\tType: " << itype << " indices: ";
562  for (const auto& flaggedelement : _splayedblock[itype]) {
563  splayout << flaggedelement->index() << ' ';
564  }
565  if (itype != _splayedblock.size() - 1)
566  splayout << std::endl;
567  }
568  LOGVERB("PFEGammaAlgo") << splayout.str();
569 #endif
570 
571  // precleaning of the ECAL clusters with respect to primary KF tracks
572  // we don't allow clusters in super clusters to be locked out this way
573  removeOrLinkECALClustersToKFTracks();
574 
575  initializeProtoCands(refinableObjects);
576  LOGDRESSED("PFEGammaAlgo") << "Initialized " << refinableObjects.size() << " proto-EGamma objects" << std::endl;
577  dumpCurrentRefinableObjects();
578 
579  //
580  // now we start the refining steps
581  //
582  //
583 
584  // --- Primary Linking Step ---
585  // since this is particle flow and we try to work from the pixels out
586  // we start by linking the tracks together and finding the ECAL clusters
587  for (auto& RO : refinableObjects) {
588  // find the KF tracks associated to GSF primary tracks
589  linkRefinableObjectGSFTracksToKFs(RO);
590  // do the same for HCAL clusters associated to the GSF
591  linkRefinableObjectPrimaryGSFTrackToHCAL(RO);
592  // link secondary KF tracks associated to primary KF tracks
593  linkRefinableObjectPrimaryKFsToSecondaryKFs(RO);
594  // pick up clusters that are linked to the GSF primary
595  linkRefinableObjectPrimaryGSFTrackToECAL(RO);
596  // link associated KF to ECAL (ECAL part grabs PS clusters too if able)
597  linkRefinableObjectKFTracksToECAL(RO);
598  // now finally look for clusters associated to brem tangents
599  linkRefinableObjectBremTangentsToECAL(RO);
600  }
601 
602  LOGDRESSED("PFEGammaAlgo") << "Dumping after GSF and KF Track (Primary) Linking : " << std::endl;
603  dumpCurrentRefinableObjects();
604 
605  // merge objects after primary linking
606  mergeROsByAnyLink(refinableObjects);
607 
608  LOGDRESSED("PFEGammaAlgo") << "Dumping after first merging operation : " << std::endl;
609  dumpCurrentRefinableObjects();
610 
611  // --- Secondary Linking Step ---
612  // after this we go through the ECAL clusters on the remaining tracks
613  // and try to link those in...
614  for (auto& RO : refinableObjects) {
615  // look for conversion legs
616  linkRefinableObjectECALToSingleLegConv(RO);
617  dumpCurrentRefinableObjects();
618  // look for tracks that complement conversion legs
619  linkRefinableObjectConvSecondaryKFsToSecondaryKFs(RO);
620  // look again for ECAL clusters (this time with an e/p cut)
621  linkRefinableObjectSecondaryKFsToECAL(RO);
622  }
623 
624  LOGDRESSED("PFEGammaAlgo") << "Dumping after ECAL to Track (Secondary) Linking : " << std::endl;
625  dumpCurrentRefinableObjects();
626 
627  // merge objects after primary linking
628  mergeROsByAnyLink(refinableObjects);
629 
630  LOGDRESSED("PFEGammaAlgo") << "There are " << refinableObjects.size() << " after the 2nd merging step." << std::endl;
631  dumpCurrentRefinableObjects();
632 
633  // -- unlinking and proto-object vetos, final sorting
634  for (auto& RO : refinableObjects) {
635  // remove secondary KFs (and possibly ECALs) matched to HCAL clusters
636  unlinkRefinableObjectKFandECALMatchedToHCAL(RO, false, false);
637  // remove secondary KFs and ECALs linked to them that have bad E/p_in
638  // and spoil the resolution
639  unlinkRefinableObjectKFandECALWithBadEoverP(RO);
640  // put things back in order after partitioning
641  std::sort(RO.ecalclusters.begin(), RO.ecalclusters.end(), [](auto const& a, auto const& b) {
642  return (a->clusterRef()->correctedEnergy() > b->clusterRef()->correctedEnergy());
643  });
644  setROElectronCluster(RO);
645  }
646 
647  LOGDRESSED("PFEGammaAlgo") << "There are " << refinableObjects.size() << " after the unlinking and vetos step."
648  << std::endl;
649  dumpCurrentRefinableObjects();
650 
651  // fill the PF candidates and then build the refined SC
652  return fillPFCandidates(refinableObjects);
653 }
654 
655 void PFEGammaAlgo::initializeProtoCands(std::list<PFEGammaAlgo::ProtoEGObject>& egobjs) {
656  // step 1: build SC based proto-candidates
657  // in the future there will be an SC Et requirement made here to control
658  // block size
659  for (auto& element : _splayedblock[PFBlockElement::SC]) {
660  LOGDRESSED("PFEGammaAlgo") << "creating SC-based proto-object" << std::endl
661  << "\tSC at index: " << element->index() << " has type: " << element->type()
662  << std::endl;
663  element.setFlag(false);
664  ProtoEGObject fromSC;
665  fromSC.nBremsWithClusters = -1;
666  fromSC.firstBrem = -1;
667  fromSC.lateBrem = -1;
668  fromSC.parentBlock = _currentblock;
669  fromSC.parentSC = docast(const PFSCElement*, element.get());
670  // splay the supercluster so we can knock out used elements
671  bool sc_success = unwrapSuperCluster(fromSC.parentSC, fromSC.ecalclusters, fromSC.ecal2ps);
672  if (sc_success) {
673  /*
674  auto ins_pos = std::lower_bound(refinableObjects.begin(),
675  refinableObjects.end(),
676  fromSC,
677  [&](const ProtoEGObject& a,
678  const ProtoEGObject& b){
679  const double a_en =
680  a.parentSC->superClusterRef()->energy();
681  const double b_en =
682  b.parentSC->superClusterRef()->energy();
683  return a_en < b_en;
684  });
685  */
686  egobjs.insert(egobjs.end(), fromSC);
687  }
688  }
689  // step 2: build GSF-seed-based proto-candidates
690  reco::GsfTrackRef gsfref_forextra;
691  reco::TrackExtraRef gsftrk_extra;
692  reco::ElectronSeedRef theseedref;
693  for (auto& element : _splayedblock[PFBlockElement::GSF]) {
694  LOGDRESSED("PFEGammaAlgo") << "creating GSF-based proto-object" << std::endl
695  << "\tGSF at index: " << element->index() << " has type: " << element->type()
696  << std::endl;
697  const PFGSFElement* elementAsGSF = docast(const PFGSFElement*, element.get());
698  if (elementAsGSF->trackType(reco::PFBlockElement::T_FROM_GAMMACONV)) {
699  continue; // for now, do not allow dedicated brems to make proto-objects
700  }
701  element.setFlag(false);
702 
703  ProtoEGObject fromGSF;
704  fromGSF.nBremsWithClusters = -1;
705  fromGSF.firstBrem = -1;
706  fromGSF.lateBrem = 0;
707  gsfref_forextra = elementAsGSF->GsftrackRef();
708  gsftrk_extra = (gsfref_forextra.isAvailable() ? gsfref_forextra->extra() : reco::TrackExtraRef());
709  theseedref = (gsftrk_extra.isAvailable() ? gsftrk_extra->seedRef().castTo<reco::ElectronSeedRef>()
710  : reco::ElectronSeedRef());
711  fromGSF.electronSeed = theseedref;
712  // exception if there's no seed
713  if (fromGSF.electronSeed.isNull() || !fromGSF.electronSeed.isAvailable()) {
714  std::stringstream gsf_err;
715  elementAsGSF->Dump(gsf_err, "\t");
716  throw cms::Exception("PFEGammaAlgo::initializeProtoCands()")
717  << "Found a GSF track with no seed! This should not happen!" << std::endl
718  << gsf_err.str() << std::endl;
719  }
720  // flag this GSF element as globally used and push back the track ref
721  // into the protocand
722  element.setFlag(false);
723  fromGSF.parentBlock = _currentblock;
724  fromGSF.primaryGSFs.push_back(elementAsGSF);
725  // add the directly matched brem tangents
726  for (auto& brem : _splayedblock[PFBlockElement::BREM]) {
727  float dist =
728  _currentblock->dist(elementAsGSF->index(), brem->index(), _currentlinks, reco::PFBlock::LINKTEST_ALL);
729  if (dist == 0.001f) {
730  const PFBremElement* eAsBrem = docast(const PFBremElement*, brem.get());
731  fromGSF.brems.push_back(eAsBrem);
732  fromGSF.localMap.insert(eAsBrem, elementAsGSF);
733  brem.setFlag(false);
734  }
735  }
736  // if this track is ECAL seeded reset links or import cluster
737  // tracker (this is pixel only, right?) driven seeds just get the GSF
738  // track associated since this only branches for ECAL Driven seeds
739  if (fromGSF.electronSeed->isEcalDriven()) {
740  // step 2a: either merge with existing ProtoEG object with SC or add
741  // SC directly to this proto EG object if not present
742  LOGDRESSED("PFEGammaAlgo") << "GSF-based proto-object is ECAL driven, merging SC-cand" << std::endl;
743  LOGVERB("PFEGammaAlgo") << "ECAL Seed Ptr: " << fromGSF.electronSeed.get()
744  << " isAvailable: " << fromGSF.electronSeed.isAvailable()
745  << " isNonnull: " << fromGSF.electronSeed.isNonnull() << std::endl;
746  SeedMatchesToProtoObject sctoseedmatch(fromGSF.electronSeed);
747  auto objsbegin = egobjs.begin();
748  auto objsend = egobjs.end();
749  // this auto is a std::list<ProtoEGObject>::iterator
750  auto clusmatch = std::find_if(objsbegin, objsend, sctoseedmatch);
751  if (clusmatch != objsend) {
752  fromGSF.parentSC = clusmatch->parentSC;
753  fromGSF.ecalclusters = std::move(clusmatch->ecalclusters);
754  fromGSF.ecal2ps = std::move(clusmatch->ecal2ps);
755  egobjs.erase(clusmatch);
756  } else if (fromGSF.electronSeed.isAvailable() && fromGSF.electronSeed.isNonnull()) {
757  // link tests in the gap region can current split a gap electron
758  // HEY THIS IS A WORK AROUND FOR A KNOWN BUG IN PFBLOCKALGO
759  // MAYBE WE SHOULD FIX IT??????????????????????????????????
760  LOGDRESSED("PFEGammaAlgo") << "Encountered the known GSF-SC splitting bug "
761  << " in PFBlockAlgo! We should really fix this!" << std::endl;
762  } else { // SC was not in a earlier proto-object
763  std::stringstream gsf_err;
764  elementAsGSF->Dump(gsf_err, "\t");
765  throw cms::Exception("PFEGammaAlgo::initializeProtoCands()")
766  << "Expected SuperCluster from ECAL driven GSF seed "
767  << "was not found in the block!" << std::endl
768  << gsf_err.str() << std::endl;
769  } // supercluster in block
770  } // is ECAL driven seed?
771 
772  egobjs.insert(egobjs.end(), fromGSF);
773  } // end loop on GSF elements of block
774 }
775 
776 bool PFEGammaAlgo::unwrapSuperCluster(const PFSCElement* thesc,
777  std::vector<FlaggedPtr<const PFClusterElement>>& ecalclusters,
778  ClusterMap& ecal2ps) {
779  ecalclusters.clear();
780  ecal2ps.clear();
781  LOGVERB("PFEGammaAlgo") << "Pointer to SC element: 0x" << std::hex << thesc << std::dec << std::endl
782  << "cleared ecalclusters and ecal2ps!" << std::endl;
783  auto ecalbegin = _splayedblock[reco::PFBlockElement::ECAL].begin();
784  auto ecalend = _splayedblock[reco::PFBlockElement::ECAL].end();
785  auto hgcalbegin = _splayedblock[reco::PFBlockElement::HGCAL].begin();
786  auto hgcalend = _splayedblock[reco::PFBlockElement::HGCAL].end();
787  if (ecalbegin == ecalend && hgcalbegin == hgcalend) {
788  LOGERR("PFEGammaAlgo::unwrapSuperCluster()") << "There are no ECAL elements in a block with imported SC!"
789  << " This is a bug we should fix this!" << std::endl;
790  return false;
791  }
792  reco::SuperClusterRef scref = thesc->superClusterRef();
793  const bool is_pf_sc = thesc->fromPFSuperCluster();
794  if (!(scref.isAvailable() && scref.isNonnull())) {
795  throw cms::Exception("PFEGammaAlgo::unwrapSuperCluster()")
796  << "SuperCluster pointed to by block element is null!" << std::endl;
797  }
798  LOGDRESSED("PFEGammaAlgo") << "Got a valid super cluster ref! 0x" << std::hex << scref.get() << std::dec << std::endl;
799  const size_t nscclusters = scref->clustersSize();
800  const size_t nscpsclusters = scref->preshowerClustersSize();
801  size_t npfpsclusters = 0;
802  size_t npfclusters = 0;
803  LOGDRESSED("PFEGammaAlgo") << "Precalculated cluster multiplicities: " << nscclusters << ' ' << nscpsclusters
804  << std::endl;
805  NotCloserToOther<reco::PFBlockElement::SC, reco::PFBlockElement::ECAL> ecalClustersInSC(_currentblock, thesc);
806  NotCloserToOther<reco::PFBlockElement::SC, reco::PFBlockElement::HGCAL> hgcalClustersInSC(_currentblock, thesc);
807  auto ecalfirstnotinsc = std::partition(ecalbegin, ecalend, ecalClustersInSC);
808  auto hgcalfirstnotinsc = std::partition(hgcalbegin, hgcalend, hgcalClustersInSC);
809  //reset the begin and end iterators
810  ecalbegin = _splayedblock[reco::PFBlockElement::ECAL].begin();
811  ecalend = _splayedblock[reco::PFBlockElement::ECAL].end();
812 
813  hgcalbegin = _splayedblock[reco::PFBlockElement::HGCAL].begin();
814  hgcalend = _splayedblock[reco::PFBlockElement::HGCAL].end();
815 
816  //get list of associated clusters by det id and energy matching
817  //(only needed when using non-pf supercluster)
818  std::vector<const ClusterElement*> safePFClusters =
819  is_pf_sc ? std::vector<const ClusterElement*>()
820  : getSCAssociatedECALsSafe(scref, _splayedblock[reco::PFBlockElement::ECAL]);
821 
822  if (ecalfirstnotinsc == ecalbegin && hgcalfirstnotinsc == hgcalbegin) {
823  LOGERR("PFEGammaAlgo::unwrapSuperCluster()") << "No associated block elements to SuperCluster!"
824  << " This is a bug we should fix!" << std::endl;
825  return false;
826  }
827  npfclusters = std::distance(ecalbegin, ecalfirstnotinsc) + std::distance(hgcalbegin, hgcalfirstnotinsc);
828  // ensure we have found the correct number of PF ecal clusters in the case
829  // that this is a PF supercluster, otherwise all bets are off
830  if (is_pf_sc && nscclusters != npfclusters) {
831  std::stringstream sc_err;
832  thesc->Dump(sc_err, "\t");
833  throw cms::Exception("PFEGammaAlgo::unwrapSuperCluster()")
834  << "The number of found ecal elements (" << nscclusters << ") in block is not the same as"
835  << " the number of ecal PF clusters reported by the PFSuperCluster"
836  << " itself (" << npfclusters << ")! This should not happen!" << std::endl
837  << sc_err.str() << std::endl;
838  }
839  for (auto ecalitr = ecalbegin; ecalitr != ecalfirstnotinsc; ++ecalitr) {
840  const PFClusterElement* elemascluster = docast(const PFClusterElement*, ecalitr->get());
841 
842  // reject clusters that really shouldn't be associated to the SC
843  // (only needed when using non-pf-supercluster)
844  if (!is_pf_sc && std::find(safePFClusters.begin(), safePFClusters.end(), elemascluster) == safePFClusters.end())
845  continue;
846 
847  //add cluster
848  ecalclusters.emplace_back(elemascluster, true);
849  //mark cluster as used
850  ecalitr->setFlag(false);
851 
852  // process the ES elements
853  // auto is a pair<Iterator,bool> here, bool is false when placing fails
854  auto emplaceresult = ecal2ps.emplace(elemascluster, ClusterMap::mapped_type());
855  if (!emplaceresult.second) {
856  std::stringstream clus_err;
857  elemascluster->Dump(clus_err, "\t");
858  throw cms::Exception("PFEGammaAlgo::unwrapSuperCluster()")
859  << "List of pointers to ECAL block elements contains non-unique items!"
860  << " This is very bad!" << std::endl
861  << "cluster ptr = 0x" << std::hex << elemascluster << std::dec << std::endl
862  << clus_err.str() << std::endl;
863  }
864  ClusterMap::mapped_type& eslist = emplaceresult.first->second;
865  npfpsclusters += attachPSClusters(elemascluster, eslist);
866  } // loop over ecal elements
867 
868  for (auto hgcalitr = hgcalbegin; hgcalitr != hgcalfirstnotinsc; ++hgcalitr) {
869  const PFClusterElement* elemascluster = docast(const PFClusterElement*, hgcalitr->get());
870 
871  // reject clusters that really shouldn't be associated to the SC
872  // (only needed when using non-pf-supercluster)
873  if (!is_pf_sc && std::find(safePFClusters.begin(), safePFClusters.end(), elemascluster) == safePFClusters.end())
874  continue;
875 
876  //add cluster
877  ecalclusters.emplace_back(elemascluster, true);
878  //mark cluster as used
879  hgcalitr->setFlag(false);
880  } // loop over ecal elements
881 
882  /*
883  if( is_pf_sc && nscpsclusters != npfpsclusters) {
884  std::stringstream sc_err;
885  thesc->Dump(sc_err,"\t");
886  throw cms::Exception("PFEGammaAlgo::unwrapSuperCluster()")
887  << "The number of found PF preshower elements ("
888  << npfpsclusters << ") in block is not the same as"
889  << " the number of preshower clusters reported by the PFSuperCluster"
890  << " itself (" << nscpsclusters << ")! This should not happen!"
891  << std::endl
892  << sc_err.str() << std::endl;
893  }
894  */
895 
896  LOGDRESSED("PFEGammaAlgo") << " Unwrapped SC has " << npfclusters << " ECAL sub-clusters"
897  << " and " << npfpsclusters << " PreShower layers 1 & 2 clusters!" << std::endl;
898  return true;
899 }
900 
901 int PFEGammaAlgo::attachPSClusters(const ClusterElement* ecalclus, ClusterMap::mapped_type& eslist) {
902  if (ecalclus->clusterRef()->layer() == PFLayer::ECAL_BARREL)
903  return 0;
904  edm::Ptr<reco::PFCluster> clusptr = refToPtr(ecalclus->clusterRef());
905  EEtoPSElement ecalkey(clusptr.key(), clusptr);
906  auto assc_ps =
907  std::equal_range(eetops_.cbegin(), eetops_.cend(), ecalkey, [](const EEtoPSElement& a, const EEtoPSElement& b) {
908  return a.first < b.first;
909  });
910  for (const auto& ps1 : _splayedblock[reco::PFBlockElement::PS1]) {
911  edm::Ptr<reco::PFCluster> temp = refToPtr(ps1->clusterRef());
912  for (auto pscl = assc_ps.first; pscl != assc_ps.second; ++pscl) {
913  if (pscl->second == temp) {
914  const ClusterElement* pstemp = docast(const ClusterElement*, ps1.get());
915  eslist.emplace_back(pstemp);
916  }
917  }
918  }
919  for (const auto& ps2 : _splayedblock[reco::PFBlockElement::PS2]) {
920  edm::Ptr<reco::PFCluster> temp = refToPtr(ps2->clusterRef());
921  for (auto pscl = assc_ps.first; pscl != assc_ps.second; ++pscl) {
922  if (pscl->second == temp) {
923  const ClusterElement* pstemp = docast(const ClusterElement*, ps2.get());
924  eslist.emplace_back(pstemp);
925  }
926  }
927  }
928  return eslist.size();
929 }
930 
931 void PFEGammaAlgo::dumpCurrentRefinableObjects() const {
932 #ifdef PFLOW_DEBUG
933  edm::LogVerbatim("PFEGammaAlgo")
934  //<< "Dumping current block: " << std::endl << *_currentblock << std::endl
935  << "Dumping " << refinableObjects.size() << " refinable objects for this block: " << std::endl;
936  for (const auto& ro : refinableObjects) {
937  std::stringstream info;
938  info << "Refinable Object:" << std::endl;
939  if (ro.parentSC) {
940  info << "\tSuperCluster element attached to object:" << std::endl << '\t';
941  ro.parentSC->Dump(info, "\t");
942  info << std::endl;
943  }
944  if (ro.electronSeed.isNonnull()) {
945  info << "\tGSF element attached to object:" << std::endl;
946  ro.primaryGSFs.front()->Dump(info, "\t");
947  info << std::endl;
948  info << "firstBrem : " << ro.firstBrem << " lateBrem : " << ro.lateBrem
949  << " nBrems with cluster : " << ro.nBremsWithClusters << std::endl;
950  ;
951  if (ro.electronClusters.size() && ro.electronClusters[0]) {
952  info << "electron cluster : ";
953  ro.electronClusters[0]->Dump(info, "\t");
954  info << std::endl;
955  } else {
956  info << " no electron cluster." << std::endl;
957  }
958  }
959  if (ro.primaryKFs.size()) {
960  info << "\tPrimary KF tracks attached to object: " << std::endl;
961  for (const auto& kf : ro.primaryKFs) {
962  kf->Dump(info, "\t");
963  info << std::endl;
964  }
965  }
966  if (ro.secondaryKFs.size()) {
967  info << "\tSecondary KF tracks attached to object: " << std::endl;
968  for (const auto& kf : ro.secondaryKFs) {
969  kf->Dump(info, "\t");
970  info << std::endl;
971  }
972  }
973  if (ro.brems.size()) {
974  info << "\tBrem tangents attached to object: " << std::endl;
975  for (const auto& brem : ro.brems) {
976  brem->Dump(info, "\t");
977  info << std::endl;
978  }
979  }
980  if (ro.ecalclusters.size()) {
981  info << "\tECAL clusters attached to object: " << std::endl;
982  for (const auto& clus : ro.ecalclusters) {
983  clus->Dump(info, "\t");
984  info << std::endl;
985  if (ro.ecal2ps.find(clus) != ro.ecal2ps.end()) {
986  for (const auto& psclus : ro.ecal2ps.at(clus)) {
987  info << "\t\t Attached PS Cluster: ";
988  psclus->Dump(info, "");
989  info << std::endl;
990  }
991  }
992  }
993  }
994  edm::LogVerbatim("PFEGammaAlgo") << info.str();
995  }
996 #endif
997 }
998 
999 // look through our KF tracks in this block and match
1000 void PFEGammaAlgo::removeOrLinkECALClustersToKFTracks() {
1001  typedef std::multimap<double, unsigned> MatchedMap;
1002  typedef const reco::PFBlockElementGsfTrack* GsfTrackElementPtr;
1003  if (_splayedblock[reco::PFBlockElement::ECAL].empty() || _splayedblock[reco::PFBlockElement::TRACK].empty())
1004  return;
1005  MatchedMap matchedGSFs, matchedECALs;
1006  std::unordered_map<GsfTrackElementPtr, MatchedMap> gsf_ecal_cache;
1007  for (auto& kftrack : _splayedblock[reco::PFBlockElement::TRACK]) {
1008  matchedGSFs.clear();
1009  _currentblock->associatedElements(
1010  kftrack->index(), _currentlinks, matchedGSFs, reco::PFBlockElement::GSF, reco::PFBlock::LINKTEST_ALL);
1011  if (matchedGSFs.empty()) { // only run this if we aren't associated to GSF
1012  LesserByDistance closestTrackToECAL(_currentblock, _currentlinks, kftrack.get());
1013  auto ecalbegin = _splayedblock[reco::PFBlockElement::ECAL].begin();
1014  auto ecalend = _splayedblock[reco::PFBlockElement::ECAL].end();
1015  std::partial_sort(ecalbegin, ecalbegin + 1, ecalend, closestTrackToECAL);
1016  auto& closestECAL = _splayedblock[reco::PFBlockElement::ECAL].front();
1017  const float dist =
1018  _currentblock->dist(kftrack->index(), closestECAL->index(), _currentlinks, reco::PFBlock::LINKTEST_ALL);
1019  bool inSC = false;
1020  for (auto& sc : _splayedblock[reco::PFBlockElement::SC]) {
1021  float dist_sc =
1022  _currentblock->dist(sc->index(), closestECAL->index(), _currentlinks, reco::PFBlock::LINKTEST_ALL);
1023  if (dist_sc != -1.0f) {
1024  inSC = true;
1025  break;
1026  }
1027  }
1028 
1029  if (dist != -1.0f && closestECAL.flag()) {
1030  bool gsflinked = false;
1031  // check that this cluster is not associated to a GSF track
1032  for (const auto& gsfflag : _splayedblock[reco::PFBlockElement::GSF]) {
1033  const reco::PFBlockElementGsfTrack* elemasgsf = docast(const reco::PFBlockElementGsfTrack*, gsfflag.get());
1034  if (elemasgsf->trackType(reco::PFBlockElement::T_FROM_GAMMACONV)) {
1035  continue; // keep clusters that have a found conversion GSF near
1036  }
1037  // make sure cache exists
1038  if (!gsf_ecal_cache.count(elemasgsf)) {
1039  matchedECALs.clear();
1040  _currentblock->associatedElements(elemasgsf->index(),
1041  _currentlinks,
1042  matchedECALs,
1043  reco::PFBlockElement::ECAL,
1044  reco::PFBlock::LINKTEST_ALL);
1045  gsf_ecal_cache.emplace(elemasgsf, matchedECALs);
1046  MatchedMap().swap(matchedECALs);
1047  }
1048  const MatchedMap& ecal_matches = gsf_ecal_cache[elemasgsf];
1049  if (!ecal_matches.empty()) {
1050  if (ecal_matches.begin()->second == closestECAL->index()) {
1051  gsflinked = true;
1052  break;
1053  }
1054  }
1055  } // loop over primary GSF tracks
1056  if (!gsflinked && !inSC) {
1057  // determine if we should remove the matched cluster
1058  const reco::PFBlockElementTrack* kfEle = docast(const reco::PFBlockElementTrack*, kftrack.get());
1059  const reco::TrackRef& trackref = kfEle->trackRef();
1060 
1061  const int nexhits = trackref->hitPattern().numberOfLostHits(HitPattern::MISSING_INNER_HITS);
1062  bool fromprimaryvertex = false;
1063  for (auto vtxtks = primaryVertex_.tracks_begin(); vtxtks != primaryVertex_.tracks_end(); ++vtxtks) {
1064  if (trackref == vtxtks->castTo<reco::TrackRef>()) {
1065  fromprimaryvertex = true;
1066  break;
1067  }
1068  } // loop over tracks in primary vertex
1069  // if associated to good non-GSF matched track remove this cluster
1070  if (PFTrackAlgoTools::isGoodForEGMPrimary(trackref->algo()) && nexhits == 0 && fromprimaryvertex) {
1071  closestECAL.setFlag(false);
1072  }
1073  }
1074  } // found a good closest ECAL match
1075  } // no GSF track matched to KF
1076  } // loop over KF elements
1077 }
1078 
1079 void PFEGammaAlgo::mergeROsByAnyLink(std::list<PFEGammaAlgo::ProtoEGObject>& ROs) {
1080  if (ROs.size() < 2)
1081  return; // nothing to do with one or zero ROs
1082  bool check_for_merge = true;
1083  while (check_for_merge) {
1084  // bugfix for early termination merging loop (15 April 2014)
1085  // check all pairwise combinations in the list
1086  // if one has a merge shuffle it to the front of the list
1087  // if there are no merges left to do we can terminate
1088  for (auto it1 = ROs.begin(); it1 != ROs.end(); ++it1) {
1089  auto find_start = it1;
1090  ++find_start;
1091  auto has_merge = std::find_if(find_start, ROs.end(), std::bind(testIfROMergableByLink, _1, *it1));
1092  if (has_merge != ROs.end() && it1 != ROs.begin()) {
1093  std::swap(*(ROs.begin()), *it1);
1094  break;
1095  }
1096  } // ensure mergables are shuffled to the front
1097  ProtoEGObject& thefront = ROs.front();
1098  auto mergestart = ROs.begin();
1099  ++mergestart;
1100  auto nomerge = std::partition(mergestart, ROs.end(), std::bind(testIfROMergableByLink, _1, thefront));
1101  if (nomerge != mergestart) {
1102  LOGDRESSED("PFEGammaAlgo::mergeROsByAnyLink()")
1103  << "Found objects " << std::distance(mergestart, nomerge) << " to merge by links to the front!" << std::endl;
1104  for (auto roToMerge = mergestart; roToMerge != nomerge; ++roToMerge) {
1105  //bugfix! L.Gray 14 Jan 2016
1106  // -- check that the front is still mergeable!
1107  if (!thefront.ecalclusters.empty() && !roToMerge->ecalclusters.empty()) {
1108  if (thefront.ecalclusters.front()->clusterRef()->layer() !=
1109  roToMerge->ecalclusters.front()->clusterRef()->layer()) {
1110  LOGWARN("PFEGammaAlgo::mergeROsByAnyLink") << "Tried to merge EB and EE clusters! Skipping!";
1111  ROs.push_back(*roToMerge);
1112  continue;
1113  }
1114  }
1115  //end bugfix
1116  thefront.ecalclusters.insert(
1117  thefront.ecalclusters.end(), roToMerge->ecalclusters.begin(), roToMerge->ecalclusters.end());
1118  thefront.ecal2ps.insert(roToMerge->ecal2ps.begin(), roToMerge->ecal2ps.end());
1119  thefront.secondaryKFs.insert(
1120  thefront.secondaryKFs.end(), roToMerge->secondaryKFs.begin(), roToMerge->secondaryKFs.end());
1121 
1122  thefront.localMap.concatenate(roToMerge->localMap);
1123  // TO FIX -> use best (E_gsf - E_clustersum)/E_GSF
1124  if (!thefront.parentSC && roToMerge->parentSC) {
1125  thefront.parentSC = roToMerge->parentSC;
1126  }
1127  if (thefront.electronSeed.isNull() && roToMerge->electronSeed.isNonnull()) {
1128  thefront.electronSeed = roToMerge->electronSeed;
1129  thefront.primaryGSFs.insert(
1130  thefront.primaryGSFs.end(), roToMerge->primaryGSFs.begin(), roToMerge->primaryGSFs.end());
1131  thefront.primaryKFs.insert(
1132  thefront.primaryKFs.end(), roToMerge->primaryKFs.begin(), roToMerge->primaryKFs.end());
1133  thefront.brems.insert(thefront.brems.end(), roToMerge->brems.begin(), roToMerge->brems.end());
1134  thefront.electronClusters = roToMerge->electronClusters;
1135  thefront.nBremsWithClusters = roToMerge->nBremsWithClusters;
1136  thefront.firstBrem = roToMerge->firstBrem;
1137  thefront.lateBrem = roToMerge->lateBrem;
1138  } else if (thefront.electronSeed.isNonnull() && roToMerge->electronSeed.isNonnull()) {
1139  LOGDRESSED("PFEGammaAlgo::mergeROsByAnyLink")
1140  << "Need to implement proper merging of two gsf candidates!" << std::endl;
1141  }
1142  }
1143  ROs.erase(mergestart, nomerge);
1144  // put the merged element in the back of the cleaned list
1145  ROs.push_back(ROs.front());
1146  ROs.pop_front();
1147  } else {
1148  check_for_merge = false;
1149  }
1150  }
1151  LOGDRESSED("PFEGammaAlgo::mergeROsByAnyLink()")
1152  << "After merging by links there are: " << ROs.size() << " refinable EGamma objects!" << std::endl;
1153 }
1154 
1155 // pull in KF tracks associated to the RO but not closer to another
1156 // NB: in initializeProtoCands() we forced the GSF tracks not to be
1157 // from a conversion, but we will leave a protection here just in
1158 // case things change in the future
1159 void PFEGammaAlgo::linkRefinableObjectGSFTracksToKFs(ProtoEGObject& RO) {
1160  constexpr reco::PFBlockElement::TrackType convType = reco::PFBlockElement::T_FROM_GAMMACONV;
1161  if (_splayedblock[reco::PFBlockElement::TRACK].empty())
1162  return;
1163  auto KFbegin = _splayedblock[reco::PFBlockElement::TRACK].begin();
1164  auto KFend = _splayedblock[reco::PFBlockElement::TRACK].end();
1165  for (auto& gsfflagged : RO.primaryGSFs) {
1166  const PFGSFElement* seedtk = gsfflagged;
1167  // don't process SC-only ROs or secondary seeded ROs
1168  if (RO.electronSeed.isNull() || seedtk->trackType(convType))
1169  continue;
1170  NotCloserToOther<reco::PFBlockElement::GSF, reco::PFBlockElement::TRACK> gsfTrackToKFs(_currentblock, seedtk);
1171  // get KF tracks not closer to another and not already used
1172  auto notlinked = std::partition(KFbegin, KFend, gsfTrackToKFs);
1173  // attach tracks and set as used
1174  for (auto kft = KFbegin; kft != notlinked; ++kft) {
1175  const PFKFElement* elemaskf = docast(const PFKFElement*, kft->get());
1176  // don't care about things that aren't primaries or directly
1177  // associated secondary tracks
1178  if (isPrimaryTrack(*elemaskf, *seedtk) && !elemaskf->trackType(convType)) {
1179  kft->setFlag(false);
1180  RO.primaryKFs.push_back(elemaskf);
1181  RO.localMap.insert(seedtk, elemaskf);
1182  } else if (elemaskf->trackType(convType)) {
1183  kft->setFlag(false);
1184  RO.secondaryKFs.push_back(elemaskf);
1185  RO.localMap.insert(seedtk, elemaskf);
1186  }
1187  } // loop on closest KFs not closer to other GSFs
1188  } // loop on GSF primaries on RO
1189 }
1190 
1191 void PFEGammaAlgo::linkRefinableObjectPrimaryKFsToSecondaryKFs(ProtoEGObject& RO) {
1192  constexpr reco::PFBlockElement::TrackType convType = reco::PFBlockElement::T_FROM_GAMMACONV;
1193  if (_splayedblock[reco::PFBlockElement::TRACK].empty())
1194  return;
1195  auto KFbegin = _splayedblock[reco::PFBlockElement::TRACK].begin();
1196  auto KFend = _splayedblock[reco::PFBlockElement::TRACK].end();
1197  for (auto& primkf : RO.primaryKFs) {
1198  // don't process SC-only ROs or secondary seeded ROs
1199  if (primkf->trackType(convType)) {
1200  throw cms::Exception("PFEGammaAlgo::linkRefinableObjectPrimaryKFsToSecondaryKFs()")
1201  << "A KF track from conversion has been assigned as a primary!!" << std::endl;
1202  }
1203  NotCloserToOther<reco::PFBlockElement::TRACK, reco::PFBlockElement::TRACK, true> kfTrackToKFs(_currentblock,
1204  primkf);
1205  // get KF tracks not closer to another and not already used
1206  auto notlinked = std::partition(KFbegin, KFend, kfTrackToKFs);
1207  // attach tracks and set as used
1208  for (auto kft = KFbegin; kft != notlinked; ++kft) {
1209  const PFKFElement* elemaskf = docast(const PFKFElement*, kft->get());
1210  // don't care about things that aren't primaries or directly
1211  // associated secondary tracks
1212  if (elemaskf->trackType(convType)) {
1213  kft->setFlag(false);
1214  RO.secondaryKFs.push_back(elemaskf);
1215  RO.localMap.insert(primkf, elemaskf);
1216  }
1217  } // loop on closest KFs not closer to other KFs
1218  } // loop on KF primaries on RO
1219 }
1220 
1221 // try to associate the tracks to cluster elements which are not used
1222 void PFEGammaAlgo::linkRefinableObjectPrimaryGSFTrackToECAL(ProtoEGObject& RO) {
1223  if (_splayedblock[reco::PFBlockElement::ECAL].empty()) {
1224  RO.electronClusters.push_back(nullptr);
1225  return;
1226  }
1227  auto ECALbegin = _splayedblock[reco::PFBlockElement::ECAL].begin();
1228  auto ECALend = _splayedblock[reco::PFBlockElement::ECAL].end();
1229  for (auto& primgsf : RO.primaryGSFs) {
1230  NotCloserToOther<reco::PFBlockElement::GSF, reco::PFBlockElement::ECAL> gsfTracksToECALs(_currentblock, primgsf);
1231  // get set of matching ecals not already in SC
1232  auto notmatched_blk = std::partition(ECALbegin, ECALend, gsfTracksToECALs);
1233  notmatched_blk =
1234  std::partition(ECALbegin, notmatched_blk, [&primgsf](auto const& x) { return compatibleEoPOut(*x, *primgsf); });
1235  // get set of matching ecals already in the RO
1236  auto notmatched_sc = std::partition(RO.ecalclusters.begin(), RO.ecalclusters.end(), gsfTracksToECALs);
1237  notmatched_sc = std::partition(
1238  RO.ecalclusters.begin(), notmatched_sc, [&primgsf](auto const& x) { return compatibleEoPOut(*x, *primgsf); });
1239  // look inside the SC for the ECAL cluster
1240  for (auto ecal = RO.ecalclusters.begin(); ecal != notmatched_sc; ++ecal) {
1241  const PFClusterElement* elemascluster = docast(const PFClusterElement*, ecal->get());
1242  FlaggedPtr<const PFClusterElement> temp(elemascluster, true);
1243  LOGDRESSED("PFEGammaAlgo::linkGSFTracktoECAL()") << "Found a cluster already in RO by GSF extrapolation"
1244  << " at ECAL surface!" << std::endl
1245  << *elemascluster << std::endl;
1246 
1247  RO.localMap.insert(primgsf, temp.get());
1248  }
1249  // look outside the SC for the ecal cluster
1250  for (auto ecal = ECALbegin; ecal != notmatched_blk; ++ecal) {
1251  const PFClusterElement* elemascluster = docast(const PFClusterElement*, ecal->get());
1252  LOGDRESSED("PFEGammaAlgo::linkGSFTracktoECAL()") << "Found a cluster not already in RO by GSF extrapolation"
1253  << " at ECAL surface!" << std::endl
1254  << *elemascluster << std::endl;
1255  if (addPFClusterToROSafe(elemascluster, RO)) {
1256  attachPSClusters(elemascluster, RO.ecal2ps[elemascluster]);
1257  RO.localMap.insert(primgsf, elemascluster);
1258  ecal->setFlag(false);
1259  }
1260  }
1261  }
1262 }
1263 
1264 // try to associate the tracks to cluster elements which are not used
1265 void PFEGammaAlgo::linkRefinableObjectPrimaryGSFTrackToHCAL(ProtoEGObject& RO) {
1266  if (_splayedblock[reco::PFBlockElement::HCAL].empty())
1267  return;
1268  auto HCALbegin = _splayedblock[reco::PFBlockElement::HCAL].begin();
1269  auto HCALend = _splayedblock[reco::PFBlockElement::HCAL].end();
1270  for (auto& primgsf : RO.primaryGSFs) {
1271  NotCloserToOther<reco::PFBlockElement::GSF, reco::PFBlockElement::HCAL> gsfTracksToHCALs(_currentblock, primgsf);
1272  auto notmatched = std::partition(HCALbegin, HCALend, gsfTracksToHCALs);
1273  for (auto hcal = HCALbegin; hcal != notmatched; ++hcal) {
1274  const PFClusterElement* elemascluster = docast(const PFClusterElement*, hcal->get());
1275  FlaggedPtr<const PFClusterElement> temp(elemascluster, true);
1276  LOGDRESSED("PFEGammaAlgo::linkGSFTracktoECAL()")
1277  << "Found an HCAL cluster associated to GSF extrapolation" << std::endl;
1278  RO.hcalClusters.push_back(temp.get());
1279  RO.localMap.insert(primgsf, temp.get());
1280  hcal->setFlag(false);
1281  }
1282  }
1283 }
1284 
1285 // try to associate the tracks to cluster elements which are not used
1286 void PFEGammaAlgo::linkRefinableObjectKFTracksToECAL(ProtoEGObject& RO) {
1287  if (_splayedblock[reco::PFBlockElement::ECAL].empty())
1288  return;
1289  for (auto& primkf : RO.primaryKFs)
1290  linkKFTrackToECAL(primkf, RO);
1291  for (auto& secdkf : RO.secondaryKFs)
1292  linkKFTrackToECAL(secdkf, RO);
1293 }
1294 
1295 void PFEGammaAlgo::linkKFTrackToECAL(PFKFElement const* kfflagged, ProtoEGObject& RO) {
1296  std::vector<FlaggedPtr<const PFClusterElement>>& currentECAL = RO.ecalclusters;
1297  auto ECALbegin = _splayedblock[reco::PFBlockElement::ECAL].begin();
1298  auto ECALend = _splayedblock[reco::PFBlockElement::ECAL].end();
1299  NotCloserToOther<reco::PFBlockElement::TRACK, reco::PFBlockElement::ECAL> kfTrackToECALs(_currentblock, kfflagged);
1300  NotCloserToOther<reco::PFBlockElement::GSF, reco::PFBlockElement::ECAL> kfTrackGSFToECALs(_currentblock, kfflagged);
1301  //get the ECAL elements not used and not closer to another KF
1302  auto notmatched_sc = std::partition(currentECAL.begin(), currentECAL.end(), kfTrackToECALs);
1303  //get subset ECAL elements not used or closer to another GSF of any type
1304  notmatched_sc = std::partition(currentECAL.begin(), notmatched_sc, kfTrackGSFToECALs);
1305  for (auto ecalitr = currentECAL.begin(); ecalitr != notmatched_sc; ++ecalitr) {
1306  const PFClusterElement* elemascluster = docast(const PFClusterElement*, ecalitr->get());
1307  FlaggedPtr<const PFClusterElement> flaggedclus(elemascluster, true);
1308 
1309  LOGDRESSED("PFEGammaAlgo::linkKFTracktoECAL()") << "Found a cluster already in RO by KF extrapolation"
1310  << " at ECAL surface!" << std::endl
1311  << *elemascluster << std::endl;
1312  RO.localMap.insert(elemascluster, kfflagged);
1313  }
1314  //get the ECAL elements not used and not closer to another KF
1315  auto notmatched_blk = std::partition(ECALbegin, ECALend, kfTrackToECALs);
1316  //get subset ECAL elements not used or closer to another GSF of any type
1317  notmatched_blk = std::partition(ECALbegin, notmatched_blk, kfTrackGSFToECALs);
1318  for (auto ecalitr = ECALbegin; ecalitr != notmatched_blk; ++ecalitr) {
1319  const PFClusterElement* elemascluster = docast(const PFClusterElement*, ecalitr->get());
1320  if (addPFClusterToROSafe(elemascluster, RO)) {
1321  attachPSClusters(elemascluster, RO.ecal2ps[elemascluster]);
1322  ecalitr->setFlag(false);
1323 
1324  LOGDRESSED("PFEGammaAlgo::linkKFTracktoECAL()") << "Found a cluster not in RO by KF extrapolation"
1325  << " at ECAL surface!" << std::endl
1326  << *elemascluster << std::endl;
1327  RO.localMap.insert(elemascluster, kfflagged);
1328  }
1329  }
1330 }
1331 
1332 void PFEGammaAlgo::linkRefinableObjectBremTangentsToECAL(ProtoEGObject& RO) {
1333  if (RO.brems.empty())
1334  return;
1335  int FirstBrem = -1;
1336  int TrajPos = -1;
1337  int lastBremTrajPos = -1;
1338  for (auto& brem : RO.brems) {
1339  bool has_clusters = false;
1340  TrajPos = (brem->indTrajPoint()) - 2;
1341  auto ECALbegin = _splayedblock[reco::PFBlockElement::ECAL].begin();
1342  auto ECALend = _splayedblock[reco::PFBlockElement::ECAL].end();
1343  NotCloserToOther<reco::PFBlockElement::BREM, reco::PFBlockElement::ECAL> BremToECALs(_currentblock, brem);
1344  // check for late brem using clusters already in the SC
1345  auto RSCBegin = RO.ecalclusters.begin();
1346  auto RSCEnd = RO.ecalclusters.end();
1347  auto notmatched_rsc = std::partition(RSCBegin, RSCEnd, BremToECALs);
1348  for (auto ecal = RSCBegin; ecal != notmatched_rsc; ++ecal) {
1349  float deta = std::abs((*ecal)->clusterRef()->positionREP().eta() - brem->positionAtECALEntrance().eta());
1350  if (deta < 0.015) {
1351  has_clusters = true;
1352  if (lastBremTrajPos == -1 || lastBremTrajPos < TrajPos) {
1353  lastBremTrajPos = TrajPos;
1354  }
1355  if (FirstBrem == -1 || TrajPos < FirstBrem) { // set brem information
1356  FirstBrem = TrajPos;
1357  RO.firstBrem = TrajPos;
1358  }
1359  LOGDRESSED("PFEGammaAlgo::linkBremToECAL()") << "Found a cluster already in SC linked to brem extrapolation"
1360  << " at ECAL surface!" << std::endl;
1361  RO.localMap.insert(ecal->get(), brem);
1362  }
1363  }
1364  // grab new clusters from the block (ensured to not be late brem)
1365  auto notmatched_block = std::partition(ECALbegin, ECALend, BremToECALs);
1366  for (auto ecal = ECALbegin; ecal != notmatched_block; ++ecal) {
1367  float deta = std::abs((*ecal)->clusterRef()->positionREP().eta() - brem->positionAtECALEntrance().eta());
1368  if (deta < 0.015) {
1369  has_clusters = true;
1370  if (lastBremTrajPos == -1 || lastBremTrajPos < TrajPos) {
1371  lastBremTrajPos = TrajPos;
1372  }
1373  if (FirstBrem == -1 || TrajPos < FirstBrem) { // set brem information
1374 
1375  FirstBrem = TrajPos;
1376  RO.firstBrem = TrajPos;
1377  }
1378  const PFClusterElement* elemasclus = docast(const PFClusterElement*, ecal->get());
1379  if (addPFClusterToROSafe(elemasclus, RO)) {
1380  attachPSClusters(elemasclus, RO.ecal2ps[elemasclus]);
1381 
1382  RO.localMap.insert(ecal->get(), brem);
1383  ecal->setFlag(false);
1384  LOGDRESSED("PFEGammaAlgo::linkBremToECAL()") << "Found a cluster not already associated by brem extrapolation"
1385  << " at ECAL surface!" << std::endl;
1386  }
1387  }
1388  }
1389  if (has_clusters) {
1390  if (RO.nBremsWithClusters == -1)
1391  RO.nBremsWithClusters = 0;
1392  ++RO.nBremsWithClusters;
1393  }
1394  }
1395 }
1396 
1397 void PFEGammaAlgo::linkRefinableObjectConvSecondaryKFsToSecondaryKFs(ProtoEGObject& RO) {
1398  auto KFbegin = _splayedblock[reco::PFBlockElement::TRACK].begin();
1399  auto KFend = _splayedblock[reco::PFBlockElement::TRACK].end();
1400  auto BeginROskfs = RO.secondaryKFs.begin();
1401  auto EndROskfs = RO.secondaryKFs.end();
1402  auto ronotconv = std::partition(BeginROskfs, EndROskfs, [](auto const& x) { return x->trackType(ConvType); });
1403  size_t convkfs_end = std::distance(BeginROskfs, ronotconv);
1404  for (size_t idx = 0; idx < convkfs_end; ++idx) {
1405  auto const& secKFs =
1406  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
1407  NotCloserToOther<reco::PFBlockElement::TRACK, reco::PFBlockElement::TRACK, true> TracksToTracks(_currentblock,
1408  secKFs[idx]);
1409  auto notmatched = std::partition(KFbegin, KFend, TracksToTracks);
1410  notmatched = std::partition(KFbegin, notmatched, [](auto const& x) { return x->trackType(ConvType); });
1411  for (auto kf = KFbegin; kf != notmatched; ++kf) {
1412  const reco::PFBlockElementTrack* elemaskf = docast(const reco::PFBlockElementTrack*, kf->get());
1413  RO.secondaryKFs.push_back(elemaskf);
1414  RO.localMap.insert(secKFs[idx], kf->get());
1415  kf->setFlag(false);
1416  }
1417  }
1418 }
1419 
1420 void PFEGammaAlgo::linkRefinableObjectECALToSingleLegConv(ProtoEGObject& RO) {
1421  auto KFbegin = _splayedblock[reco::PFBlockElement::TRACK].begin();
1422  auto KFend = _splayedblock[reco::PFBlockElement::TRACK].end();
1423  for (auto& ecal : RO.ecalclusters) {
1424  NotCloserToOther<reco::PFBlockElement::ECAL, reco::PFBlockElement::TRACK, true> ECALToTracks(_currentblock,
1425  ecal.get());
1426  auto notmatchedkf = std::partition(KFbegin, KFend, ECALToTracks);
1427  auto notconvkf = std::partition(KFbegin, notmatchedkf, [](auto const& x) { return x->trackType(ConvType); });
1428  // add identified KF conversion tracks
1429  for (auto kf = KFbegin; kf != notconvkf; ++kf) {
1430  const reco::PFBlockElementTrack* elemaskf = docast(const reco::PFBlockElementTrack*, kf->get());
1431  RO.secondaryKFs.push_back(elemaskf);
1432  RO.localMap.insert(ecal.get(), elemaskf);
1433  kf->setFlag(false);
1434  }
1435  // go through non-conv-identified kfs and check MVA to add conversions
1436  for (auto kf = notconvkf; kf != notmatchedkf; ++kf) {
1437  float mvaval = evaluateSingleLegMVA(_currentblock, primaryVertex_, (*kf)->index());
1438  if (mvaval > cfg_.mvaConvCut) {
1439  const reco::PFBlockElementTrack* elemaskf = docast(const reco::PFBlockElementTrack*, kf->get());
1440  RO.secondaryKFs.push_back(elemaskf);
1441  RO.localMap.insert(ecal.get(), elemaskf);
1442  kf->setFlag(false);
1443 
1444  RO.singleLegConversionMvaMap.emplace(elemaskf, mvaval);
1445  }
1446  }
1447  }
1448 }
1449 
1450 void PFEGammaAlgo::linkRefinableObjectSecondaryKFsToECAL(ProtoEGObject& RO) {
1451  auto ECALbegin = _splayedblock[reco::PFBlockElement::ECAL].begin();
1452  auto ECALend = _splayedblock[reco::PFBlockElement::ECAL].end();
1453  for (auto& skf : RO.secondaryKFs) {
1454  NotCloserToOther<reco::PFBlockElement::TRACK, reco::PFBlockElement::ECAL, false> TracksToECALwithCut(
1455  _currentblock, skf, 1.5f);
1456  auto notmatched = std::partition(ECALbegin, ECALend, TracksToECALwithCut);
1457  for (auto ecal = ECALbegin; ecal != notmatched; ++ecal) {
1458  const reco::PFBlockElementCluster* elemascluster = docast(const reco::PFBlockElementCluster*, ecal->get());
1459  if (addPFClusterToROSafe(elemascluster, RO)) {
1460  attachPSClusters(elemascluster, RO.ecal2ps[elemascluster]);
1461  RO.localMap.insert(skf, elemascluster);
1462  ecal->setFlag(false);
1463  }
1464  }
1465  }
1466 }
1467 
1468 PFEGammaAlgo::EgammaObjects PFEGammaAlgo::fillPFCandidates(const std::list<PFEGammaAlgo::ProtoEGObject>& ROs) {
1469  EgammaObjects output;
1470 
1471  // reserve output collections
1472  output.candidates.reserve(ROs.size());
1473  output.candidateExtras.reserve(ROs.size());
1474  output.refinedSuperClusters.reserve(ROs.size());
1475 
1476  for (auto& RO : ROs) {
1477  if (RO.ecalclusters.empty() && !cfg_.produceEGCandsWithNoSuperCluster)
1478  continue;
1479 
1480  reco::PFCandidate cand;
1481  reco::PFCandidateEGammaExtra xtra;
1482  if (!RO.primaryGSFs.empty() || !RO.primaryKFs.empty()) {
1483  cand.setPdgId(-11); // anything with a primary track is an electron
1484  } else {
1485  cand.setPdgId(22); // anything with no primary track is a photon
1486  }
1487  if (!RO.primaryKFs.empty()) {
1488  cand.setCharge(RO.primaryKFs[0]->trackRef()->charge());
1489  xtra.setKfTrackRef(RO.primaryKFs[0]->trackRef());
1490  cand.setTrackRef(RO.primaryKFs[0]->trackRef());
1491  cand.addElementInBlock(_currentblock, RO.primaryKFs[0]->index());
1492  }
1493  if (!RO.primaryGSFs.empty()) {
1494  cand.setCharge(RO.primaryGSFs[0]->GsftrackRef()->chargeMode());
1495  xtra.setGsfTrackRef(RO.primaryGSFs[0]->GsftrackRef());
1496  cand.setGsfTrackRef(RO.primaryGSFs[0]->GsftrackRef());
1497  cand.addElementInBlock(_currentblock, RO.primaryGSFs[0]->index());
1498  }
1499  if (RO.parentSC) {
1500  xtra.setSuperClusterPFECALRef(RO.parentSC->superClusterRef());
1501  // we'll set to the refined supercluster back up in the producer
1502  cand.setSuperClusterRef(RO.parentSC->superClusterRef());
1503  xtra.setSuperClusterRef(RO.parentSC->superClusterRef());
1504  cand.addElementInBlock(_currentblock, RO.parentSC->index());
1505  }
1506  // add brems
1507  for (const auto& brem : RO.brems) {
1508  cand.addElementInBlock(_currentblock, brem->index());
1509  }
1510  // add clusters and ps clusters
1511  for (const auto& ecal : RO.ecalclusters) {
1512  const PFClusterElement* clus = ecal.get();
1513  cand.addElementInBlock(_currentblock, clus->index());
1514  if (RO.ecal2ps.count(clus)) {
1515  for (auto& psclus : RO.ecal2ps.at(clus)) {
1516  cand.addElementInBlock(_currentblock, psclus->index());
1517  }
1518  }
1519  }
1520  // add secondary tracks
1521  for (const auto& kf : RO.secondaryKFs) {
1522  cand.addElementInBlock(_currentblock, kf->index());
1523  const reco::ConversionRefVector& convrefs = kf->convRefs();
1524  bool no_conv_ref = true;
1525  for (const auto& convref : convrefs) {
1526  if (convref.isNonnull() && convref.isAvailable()) {
1527  xtra.addConversionRef(convref);
1528  no_conv_ref = false;
1529  }
1530  }
1531  if (no_conv_ref) {
1532  //single leg conversions
1533 
1534  //look for stored mva value in map or else recompute
1535  const auto& mvavalmapped = RO.singleLegConversionMvaMap.find(kf);
1536  //FIXME: Abuse single mva value to store both provenance and single leg mva score
1537  //by storing 3.0 + mvaval
1538  float mvaval = (mvavalmapped != RO.singleLegConversionMvaMap.end()
1539  ? mvavalmapped->second
1540  : 3.0 + evaluateSingleLegMVA(_currentblock, primaryVertex_, kf->index()));
1541 
1542  xtra.addSingleLegConvTrackRefMva(std::make_pair(kf->trackRef(), mvaval));
1543  }
1544  }
1545 
1546  // build the refined supercluster from those clusters left in the cand
1547  output.refinedSuperClusters.push_back(buildRefinedSuperCluster(RO));
1548 
1549  //*TODO* cluster time is not reliable at the moment, so only use track timing
1550  float trkTime = 0, trkTimeErr = -1;
1551  if (!RO.primaryGSFs.empty() && RO.primaryGSFs[0]->isTimeValid()) {
1552  trkTime = RO.primaryGSFs[0]->time();
1553  trkTimeErr = RO.primaryGSFs[0]->timeError();
1554  } else if (!RO.primaryKFs.empty() && RO.primaryKFs[0]->isTimeValid()) {
1555  trkTime = RO.primaryKFs[0]->time();
1556  trkTimeErr = RO.primaryKFs[0]->timeError();
1557  }
1558  if (trkTimeErr >= 0) {
1559  cand.setTime(trkTime, trkTimeErr);
1560  }
1561 
1562  const reco::SuperCluster& the_sc = output.refinedSuperClusters.back();
1563  // with the refined SC in hand we build a naive candidate p4
1564  // and set the candidate ECAL position to either the barycenter of the
1565  // supercluster (if super-cluster present) or the seed of the
1566  // new SC generated by the EGAlgo
1567  const double scE = the_sc.energy();
1568  if (scE != 0.0) {
1569  const math::XYZPoint& seedPos = the_sc.seed()->position();
1570  math::XYZVector egDir = the_sc.position() - primaryVertex_.position();
1571  egDir = egDir.Unit();
1572  cand.setP4(math::XYZTLorentzVector(scE * egDir.x(), scE * egDir.y(), scE * egDir.z(), scE));
1573  math::XYZPointF ecalPOS_f(seedPos.x(), seedPos.y(), seedPos.z());
1574  cand.setPositionAtECALEntrance(ecalPOS_f);
1575  cand.setEcalEnergy(the_sc.rawEnergy(), the_sc.energy());
1576  } else if (cfg_.produceEGCandsWithNoSuperCluster && !RO.primaryGSFs.empty()) {
1577  const PFGSFElement* gsf = RO.primaryGSFs[0];
1578  const reco::GsfTrackRef& gref = gsf->GsftrackRef();
1579  math::XYZTLorentzVector p4(gref->pxMode(), gref->pyMode(), gref->pzMode(), gref->pMode());
1580  cand.setP4(p4);
1581  cand.setPositionAtECALEntrance(gsf->positionAtECALEntrance());
1582  } else if (cfg_.produceEGCandsWithNoSuperCluster && !RO.primaryKFs.empty()) {
1583  const PFKFElement* kf = RO.primaryKFs[0];
1584  reco::TrackRef kref = RO.primaryKFs[0]->trackRef();
1585  math::XYZTLorentzVector p4(kref->px(), kref->py(), kref->pz(), kref->p());
1586  cand.setP4(p4);
1587  cand.setPositionAtECALEntrance(kf->positionAtECALEntrance());
1588  }
1589  const float eleMVAValue = calculateEleMVA(RO, xtra);
1590  fillExtraInfo(RO, xtra);
1591  //std::cout << "PFEG eleMVA: " << eleMVAValue << std::endl;
1592  xtra.setMVA(eleMVAValue);
1593  cand.set_mva_e_pi(eleMVAValue);
1594  output.candidates.push_back(cand);
1595  output.candidateExtras.push_back(xtra);
1596  }
1597 
1598  return output;
1599 }
1600 
1601 float PFEGammaAlgo::calculateEleMVA(const PFEGammaAlgo::ProtoEGObject& ro, reco::PFCandidateEGammaExtra& xtra) const {
1602  if (ro.primaryGSFs.empty()) {
1603  return -2.0f;
1604  }
1605  const PFGSFElement* gsfElement = ro.primaryGSFs.front();
1606  const PFKFElement* kfElement = nullptr;
1607  if (!ro.primaryKFs.empty()) {
1608  kfElement = ro.primaryKFs.front();
1609  }
1610  auto const& refGsf = gsfElement->GsftrackRef();
1611  reco::TrackRef refKf;
1612  constexpr float mEl = 0.000511;
1613  const double eInGsf = std::hypot(refGsf->pMode(), mEl);
1614  double dEtGsfEcal = 1e6;
1615  double sigmaEtaEta = 1e-14;
1616  const double eneHcalGsf =
1617  std::accumulate(ro.hcalClusters.begin(), ro.hcalClusters.end(), 0.0, [](const double a, auto const& b) {
1618  return a + b->clusterRef()->energy();
1619  });
1620  if (!ro.primaryKFs.empty()) {
1621  refKf = ro.primaryKFs.front()->trackRef();
1622  }
1623  const double eOutGsf = gsfElement->Pout().t();
1624  const double etaOutGsf = gsfElement->positionAtECALEntrance().eta();
1625  double firstEcalGsfEnergy{0.0};
1626  double otherEcalGsfEnergy{0.0};
1627  double ecalBremEnergy{0.0};
1628  //shower shape of cluster closest to gsf track
1629  std::vector<const reco::PFCluster*> gsfCluster;
1630  for (const auto& ecal : ro.ecalclusters) {
1631  const double cenergy = ecal->clusterRef()->correctedEnergy();
1632  bool hasgsf = ro.localMap.contains(gsfElement, ecal.get());
1633  bool haskf = ro.localMap.contains(kfElement, ecal.get());
1634  bool hasbrem = false;
1635  for (const auto& brem : ro.brems) {
1636  if (ro.localMap.contains(brem, ecal.get())) {
1637  hasbrem = true;
1638  }
1639  }
1640  if (hasbrem && ecal.get() != ro.electronClusters[0]) {
1641  ecalBremEnergy += cenergy;
1642  }
1643  if (!hasbrem && ecal.get() != ro.electronClusters[0]) {
1644  if (hasgsf)
1645  otherEcalGsfEnergy += cenergy;
1646  if (haskf)
1647  ecalBremEnergy += cenergy; // from conv. brem!
1648  if (!(hasgsf || haskf))
1649  otherEcalGsfEnergy += cenergy; // stuff from SC
1650  }
1651  }
1652 
1653  if (ro.electronClusters[0]) {
1654  reco::PFClusterRef cref = ro.electronClusters[0]->clusterRef();
1655  xtra.setGsfElectronClusterRef(_currentblock, *(ro.electronClusters[0]));
1656  firstEcalGsfEnergy = cref->correctedEnergy();
1657  dEtGsfEcal = cref->positionREP().eta() - etaOutGsf;
1658  gsfCluster.push_back(cref.get());
1659  PFClusterWidthAlgo pfwidth(gsfCluster);
1660  sigmaEtaEta = pfwidth.pflowSigmaEtaEta();
1661  }
1662 
1663  // brem sequence information
1664  float firstBrem{-1.0f};
1665  float earlyBrem{-1.0f};
1666  float lateBrem{-1.0f};
1667  if (ro.nBremsWithClusters > 0) {
1668  firstBrem = ro.firstBrem;
1669  earlyBrem = ro.firstBrem < 4 ? 1.0f : 0.0f;
1670  lateBrem = ro.lateBrem == 1 ? 1.0f : 0.0f;
1671  }
1672  xtra.setEarlyBrem(earlyBrem);
1673  xtra.setLateBrem(lateBrem);
1674  if (firstEcalGsfEnergy > 0.0) {
1675  if (refGsf.isNonnull()) {
1676  xtra.setGsfTrackPout(gsfElement->Pout());
1677  // normalization observables
1678  const float ptGsf = refGsf->ptMode();
1679  const float etaGsf = refGsf->etaMode();
1680  // tracking observables
1681  const double ptModeErrorGsf = refGsf->ptModeError();
1682  float ptModeErrOverPtGsf = (ptModeErrorGsf > 0. ? ptModeErrorGsf / ptGsf : 1.0);
1683  float chi2Gsf = refGsf->normalizedChi2();
1684  float dPtOverPtGsf = (ptGsf - gsfElement->Pout().pt()) / ptGsf;
1685  // kalman filter vars
1686  float nHitKf = refKf.isNonnull() ? refKf->hitPattern().trackerLayersWithMeasurement() : 0;
1687  float chi2Kf = refKf.isNonnull() ? refKf->normalizedChi2() : -0.01;
1688 
1689  //tracker + calorimetry observables
1690  float eTotPinMode = (firstEcalGsfEnergy + otherEcalGsfEnergy + ecalBremEnergy) / eInGsf;
1691  float eGsfPoutMode = firstEcalGsfEnergy / eOutGsf;
1692  float eTotBremPinPoutMode = (ecalBremEnergy + otherEcalGsfEnergy) / (eInGsf - eOutGsf);
1693  float dEtaGsfEcalClust = std::abs(dEtGsfEcal);
1694  float logSigmaEtaEta = std::log(sigmaEtaEta);
1695  float hOverHe = eneHcalGsf / (eneHcalGsf + firstEcalGsfEnergy);
1696 
1697  xtra.setDeltaEta(dEtaGsfEcalClust);
1698  xtra.setSigmaEtaEta(sigmaEtaEta);
1699  xtra.setHadEnergy(eneHcalGsf);
1700 
1701  // Apply bounds to variables and calculate MVA
1702  dPtOverPtGsf = std::clamp(dPtOverPtGsf, -0.2f, 1.0f);
1703  ptModeErrOverPtGsf = std::min(ptModeErrOverPtGsf, 0.3f);
1704  chi2Gsf = std::min(chi2Gsf, 10.0f);
1705  chi2Kf = std::min(chi2Kf, 10.0f);
1706  eTotPinMode = std::clamp(eTotPinMode, 0.0f, 5.0f);
1707  eGsfPoutMode = std::clamp(eGsfPoutMode, 0.0f, 5.0f);
1708  eTotBremPinPoutMode = std::clamp(eTotBremPinPoutMode, 0.0f, 5.0f);
1709  dEtaGsfEcalClust = std::min(dEtaGsfEcalClust, 0.1f);
1710  logSigmaEtaEta = std::max(logSigmaEtaEta, -14.0f);
1711 
1712  // not used for moment, weird behavior of variable
1713  //float dPtOverPtKf = refKf.isNonnull() ? (refKf->pt() - refKf->outerPt())/refKf->pt() : -0.01;
1714  //dPtOverPtKf = std::clamp(dPtOverPtKf,-0.2f, 1.0f);
1715 
1716  /*
1717  * To be used for debugging:
1718  * pretty-print the PFEgamma electron MVA input variables
1719  *
1720  * std::cout << " **** PFEG BDT observables ****" << endl;
1721  * std::cout << " < Normalization > " << endl;
1722  * std::cout << " ptGsf " << ptGsf << " Pin " << eInGsf
1723  * << " Pout " << eOutGsf << " etaGsf " << etaGsf << endl;
1724  * std::cout << " < PureTracking > " << endl;
1725  * std::cout << " ptModeErrOverPtGsf " << ptModeErrOverPtGsf
1726  * << " dPtOverPtGsf " << dPtOverPtGsf
1727  * << " chi2Gsf " << chi2Gsf
1728  * << " nhit_gsf " << nhit_gsf
1729  * << " dPtOverPtKf " << dPtOverPtKf
1730  * << " chi2Kf " << chi2Kf
1731  * << " nHitKf " << nHitKf << endl;
1732  * std::cout << " < track-ecal-hcal-ps " << endl;
1733  * std::cout << " eTotPinMode " << eTotPinMode
1734  * << " eGsfPoutMode " << eGsfPoutMode
1735  * << " eTotBremPinPoutMode " << eTotBremPinPoutMode
1736  * << " dEtaGsfEcalClust " << dEtaGsfEcalClust
1737  * << " logSigmaEtaEta " << logSigmaEtaEta
1738  * << " hOverHe " << hOverHe << " Hcal energy " << eneHcalGsf
1739  * << " lateBrem " << lateBrem
1740  * << " firstBrem " << firstBrem << endl;
1741  */
1742 
1743  float vars[] = {std::log(ptGsf),
1744  etaGsf,
1745  ptModeErrOverPtGsf,
1746  dPtOverPtGsf,
1747  chi2Gsf,
1748  nHitKf,
1749  chi2Kf,
1750  eTotPinMode,
1751  eGsfPoutMode,
1752  eTotBremPinPoutMode,
1753  dEtaGsfEcalClust,
1754  logSigmaEtaEta,
1755  hOverHe,
1756  lateBrem,
1757  firstBrem};
1758 
1759  return gbrForests_.ele_->GetAdaBoostClassifier(vars);
1760  }
1761  }
1762  return -2.0f;
1763 }
1764 
1765 void PFEGammaAlgo::fillExtraInfo(const ProtoEGObject& RO, reco::PFCandidateEGammaExtra& xtra) {
1766  // add tracks associated to clusters that are not T_FROM_GAMMACONV
1767  // info about single-leg convs is already save, so just veto in loops
1768  auto KFbegin = _splayedblock[reco::PFBlockElement::TRACK].begin();
1769  auto KFend = _splayedblock[reco::PFBlockElement::TRACK].end();
1770  for (auto& ecal : RO.ecalclusters) {
1771  NotCloserToOther<reco::PFBlockElement::ECAL, reco::PFBlockElement::TRACK, true> ECALToTracks(_currentblock,
1772  ecal.get());
1773  auto notmatchedkf = std::partition(KFbegin, KFend, ECALToTracks);
1774  auto notconvkf = std::partition(KFbegin, notmatchedkf, [](auto const& x) { return x->trackType(ConvType); });
1775  // go through non-conv-identified kfs and check MVA to add conversions
1776  for (auto kf = notconvkf; kf != notmatchedkf; ++kf) {
1777  const reco::PFBlockElementTrack* elemaskf = docast(const reco::PFBlockElementTrack*, kf->get());
1778  xtra.addExtraNonConvTrack(_currentblock, *elemaskf);
1779  }
1780  }
1781 }
1782 
1783 // currently stolen from PFECALSuperClusterAlgo, we should
1784 // try to factor this correctly since the operation is the same in
1785 // both places...
1786 reco::SuperCluster PFEGammaAlgo::buildRefinedSuperCluster(const PFEGammaAlgo::ProtoEGObject& RO) {
1787  if (RO.ecalclusters.empty()) {
1788  return reco::SuperCluster(0.0, math::XYZPoint(0, 0, 0));
1789  }
1790 
1791  bool isEE = false;
1792  // need the vector of raw pointers for a PF width class
1793  std::vector<const reco::PFCluster*> bare_ptrs;
1794  // calculate necessary parameters and build the SC
1795  double posX(0), posY(0), posZ(0), rawSCEnergy(0), corrSCEnergy(0), corrPSEnergy(0), ps1_energy(0.0), ps2_energy(0.0);
1796  for (auto& clus : RO.ecalclusters) {
1797  double ePS1 = 0;
1798  double ePS2 = 0;
1799  isEE = PFLayer::ECAL_ENDCAP == clus->clusterRef()->layer();
1800  auto clusptr = edm::refToPtr<reco::PFClusterCollection>(clus->clusterRef());
1801  bare_ptrs.push_back(clusptr.get());
1802 
1803  const double cluseraw = clusptr->energy();
1804  double cluscalibe = clusptr->correctedEnergy();
1805  const math::XYZPoint& cluspos = clusptr->position();
1806  posX += cluseraw * cluspos.X();
1807  posY += cluseraw * cluspos.Y();
1808  posZ += cluseraw * cluspos.Z();
1809  // update EE calibrated super cluster energies
1810  if (isEE && RO.ecal2ps.count(clus.get())) {
1811  const auto& psclusters = RO.ecal2ps.at(clus.get());
1812 
1813  std::vector<reco::PFCluster const*> psClusterPointers;
1814  for (auto const& psc : psclusters) {
1815  psClusterPointers.push_back(psc->clusterRef().get());
1816  }
1817  auto calibratedEnergies = thePFEnergyCalibration_.calibrateEndcapClusterEnergies(
1818  *clusptr, psClusterPointers, channelStatus_, cfg_.applyCrackCorrections);
1819  cluscalibe = calibratedEnergies.clusterEnergy;
1820  ePS1 = calibratedEnergies.ps1Energy;
1821  ePS2 = calibratedEnergies.ps2Energy;
1822  }
1823  if (ePS1 == -1.)
1824  ePS1 = 0;
1825  if (ePS2 == -1.)
1826  ePS2 = 0;
1827 
1828  rawSCEnergy += cluseraw;
1829  corrSCEnergy += cluscalibe;
1830  ps1_energy += ePS1;
1831  ps2_energy += ePS2;
1832  corrPSEnergy += ePS1 + ePS2;
1833  }
1834  posX /= rawSCEnergy;
1835  posY /= rawSCEnergy;
1836  posZ /= rawSCEnergy;
1837 
1838  // now build the supercluster
1839  reco::SuperCluster new_sc(corrSCEnergy, math::XYZPoint(posX, posY, posZ));
1840 
1841  auto clusptr = edm::refToPtr<reco::PFClusterCollection>(RO.ecalclusters.front()->clusterRef());
1842  new_sc.setCorrectedEnergy(corrSCEnergy);
1843  new_sc.setSeed(clusptr);
1844  new_sc.setPreshowerEnergyPlane1(ps1_energy);
1845  new_sc.setPreshowerEnergyPlane2(ps2_energy);
1846  new_sc.setPreshowerEnergy(corrPSEnergy);
1847  for (const auto& clus : RO.ecalclusters) {
1848  clusptr = edm::refToPtr<reco::PFClusterCollection>(clus->clusterRef());
1849  new_sc.addCluster(clusptr);
1850  auto& hits_and_fractions = clusptr->hitsAndFractions();
1851  for (auto& hit_and_fraction : hits_and_fractions) {
1852  new_sc.addHitAndFraction(hit_and_fraction.first, hit_and_fraction.second);
1853  }
1854  // put the preshower stuff back in later
1855  if (RO.ecal2ps.count(clus.get())) {
1856  const auto& cluspsassociation = RO.ecal2ps.at(clus.get());
1857  // EE rechits should be uniquely matched to sets of pre-shower
1858  // clusters at this point, so we throw an exception if otherwise
1859  // now wrapped in EDM debug flags
1860  for (const auto& pscluselem : cluspsassociation) {
1861  edm::Ptr<reco::PFCluster> psclus = edm::refToPtr<reco::PFClusterCollection>(pscluselem->clusterRef());
1862 #ifdef PFFLOW_DEBUG
1863  auto found_pscluster =
1864  std::find(new_sc.preshowerClustersBegin(), new_sc.preshowerClustersEnd(), reco::CaloClusterPtr(psclus));
1865  if (found_pscluster == new_sc.preshowerClustersEnd()) {
1866 #endif
1867  new_sc.addPreshowerCluster(psclus);
1868 #ifdef PFFLOW_DEBUG
1869  } else {
1870  throw cms::Exception("PFECALSuperClusterAlgo::buildSuperCluster")
1871  << "Found a PS cluster matched to more than one EE cluster!" << std::endl
1872  << std::hex << psclus.get() << " == " << found_pscluster->get() << std::dec << std::endl;
1873  }
1874 #endif
1875  }
1876  }
1877  }
1878 
1879  // calculate linearly weighted cluster widths
1880  PFClusterWidthAlgo pfwidth(bare_ptrs);
1881  new_sc.setEtaWidth(pfwidth.pflowEtaWidth());
1882  new_sc.setPhiWidth(pfwidth.pflowPhiWidth());
1883 
1884  // cache the value of the raw energy
1885  new_sc.rawEnergy();
1886 
1887  return new_sc;
1888 }
1889 
1890 void PFEGammaAlgo::unlinkRefinableObjectKFandECALWithBadEoverP(ProtoEGObject& RO) {
1891  // this only means something for ROs with a primary GSF track
1892  if (RO.primaryGSFs.empty())
1893  return;
1894  // need energy sums to tell if we've added crap or not
1895  const double Pin_gsf = RO.primaryGSFs.front()->GsftrackRef()->pMode();
1896  const double gsfOuterEta = RO.primaryGSFs.front()->positionAtECALEntrance().Eta();
1897  double tot_ecal = 0.0;
1898  std::vector<double> min_brem_dists;
1899  std::vector<double> closest_brem_eta;
1900  // first get the total ecal energy (we should replace this with a cache)
1901  for (const auto& ecal : RO.ecalclusters) {
1902  tot_ecal += ecal->clusterRef()->correctedEnergy();
1903  // we also need to look at the minimum distance to brems
1904  // since energetic brems will be closer to the brem than the track
1905  double min_brem_dist = 5000.0;
1906  double eta = -999.0;
1907  for (const auto& brem : RO.brems) {
1908  const float dist = _currentblock->dist(brem->index(), ecal->index(), _currentlinks, reco::PFBlock::LINKTEST_ALL);
1909  if (dist < min_brem_dist && dist != -1.0f) {
1910  min_brem_dist = dist;
1911  eta = brem->positionAtECALEntrance().Eta();
1912  }
1913  }
1914  min_brem_dists.push_back(min_brem_dist);
1915  closest_brem_eta.push_back(eta);
1916  }
1917 
1918  // loop through the ECAL clusters and remove ECAL clusters matched to
1919  // secondary track either in *or* out of the SC if the E/pin is bad
1920  for (auto secd_kf = RO.secondaryKFs.begin(); secd_kf != RO.secondaryKFs.end(); ++secd_kf) {
1921  reco::TrackRef trkRef = (*secd_kf)->trackRef();
1922  const float secpin = (*secd_kf)->trackRef()->p();
1923  bool remove_this_kf = false;
1924  for (auto ecal = RO.ecalclusters.begin(); ecal != RO.ecalclusters.end(); ++ecal) {
1925  size_t bremidx = std::distance(RO.ecalclusters.begin(), ecal);
1926  const float minbremdist = min_brem_dists[bremidx];
1927  const double ecalenergy = (*ecal)->clusterRef()->correctedEnergy();
1928  const double Epin = ecalenergy / secpin;
1929  const double detaGsf = std::abs(gsfOuterEta - (*ecal)->clusterRef()->positionREP().Eta());
1930  const double detaBrem = std::abs(closest_brem_eta[bremidx] - (*ecal)->clusterRef()->positionREP().Eta());
1931 
1932  bool kf_matched = RO.localMap.contains(ecal->get(), *secd_kf);
1933 
1934  const float tkdist =
1935  _currentblock->dist((*secd_kf)->index(), (*ecal)->index(), _currentlinks, reco::PFBlock::LINKTEST_ALL);
1936 
1937  // do not reject this track if it is closer to a brem than the
1938  // secondary track, or if it lies in the delta-eta plane with the
1939  // gsf track or if it is in the dEta plane with the brems
1940  if (Epin > 3 && kf_matched && tkdist != -1.0f && tkdist < minbremdist && detaGsf > 0.05 && detaBrem > 0.015) {
1941  double res_with = std::abs((tot_ecal - Pin_gsf) / Pin_gsf);
1942  double res_without = std::abs((tot_ecal - ecalenergy - Pin_gsf) / Pin_gsf);
1943  if (res_without < res_with) {
1944  LOGDRESSED("PFEGammaAlgo") << " REJECTED_RES totenergy " << tot_ecal << " Pin_gsf " << Pin_gsf
1945  << " cluster to secondary " << ecalenergy << " res_with " << res_with
1946  << " res_without " << res_without << std::endl;
1947  tot_ecal -= ecalenergy;
1948  remove_this_kf = true;
1949  ecal = RO.ecalclusters.erase(ecal);
1950  if (ecal == RO.ecalclusters.end())
1951  break;
1952  }
1953  }
1954  }
1955  if (remove_this_kf) {
1956  secd_kf = RO.secondaryKFs.erase(secd_kf);
1957  if (secd_kf == RO.secondaryKFs.end())
1958  break;
1959  }
1960  }
1961 }
1962 
1963 void PFEGammaAlgo::unlinkRefinableObjectKFandECALMatchedToHCAL(ProtoEGObject& RO,
1964  bool removeFreeECAL,
1965  bool removeSCEcal) {
1966  std::vector<bool> cluster_in_sc;
1967  auto ecal_begin = RO.ecalclusters.begin();
1968  auto ecal_end = RO.ecalclusters.end();
1969  auto hcal_begin = _splayedblock[reco::PFBlockElement::HCAL].begin();
1970  auto hcal_end = _splayedblock[reco::PFBlockElement::HCAL].end();
1971  for (auto secd_kf = RO.secondaryKFs.begin(); secd_kf != RO.secondaryKFs.end(); ++secd_kf) {
1972  bool remove_this_kf = false;
1973  NotCloserToOther<reco::PFBlockElement::TRACK, reco::PFBlockElement::HCAL> tracksToHCALs(_currentblock, *secd_kf);
1974  reco::TrackRef trkRef = (*secd_kf)->trackRef();
1975 
1976  bool goodTrack = PFTrackAlgoTools::isGoodForEGM(trkRef->algo());
1977  const float secpin = trkRef->p();
1978 
1979  for (auto ecal = ecal_begin; ecal != ecal_end; ++ecal) {
1980  const double ecalenergy = (*ecal)->clusterRef()->correctedEnergy();
1981  // first check if the cluster is in the SC (use dist calc for fastness)
1982  const size_t clus_idx = std::distance(ecal_begin, ecal);
1983  if (cluster_in_sc.size() < clus_idx + 1) {
1984  float dist = -1.0f;
1985  if (RO.parentSC) {
1986  dist = _currentblock->dist((*secd_kf)->index(), (*ecal)->index(), _currentlinks, reco::PFBlock::LINKTEST_ALL);
1987  }
1988  cluster_in_sc.push_back(dist != -1.0f);
1989  }
1990 
1991  // if we've found a secondary KF that matches this ecal cluster
1992  // now we see if it is matched to HCAL
1993  // if it is matched to an HCAL cluster we take different
1994  // actions if the cluster was in an SC or not
1995  if (RO.localMap.contains(ecal->get(), *secd_kf)) {
1996  auto hcal_matched = std::partition(hcal_begin, hcal_end, tracksToHCALs);
1997  for (auto hcalclus = hcal_begin; hcalclus != hcal_matched; ++hcalclus) {
1998  const reco::PFBlockElementCluster* clusthcal =
1999  dynamic_cast<const reco::PFBlockElementCluster*>(hcalclus->get());
2000  const double hcalenergy = clusthcal->clusterRef()->energy();
2001  const double hpluse = ecalenergy + hcalenergy;
2002  const bool isHoHE = ((hcalenergy / hpluse) > 0.1 && goodTrack);
2003  const bool isHoE = (hcalenergy > ecalenergy);
2004  const bool isPoHE = (secpin > hpluse);
2005  if (cluster_in_sc[clus_idx]) {
2006  if (isHoE || isPoHE) {
2007  LOGDRESSED("PFEGammaAlgo") << "REJECTED TRACK FOR H/E or P/(H+E), CLUSTER IN SC"
2008  << " H/H+E " << (hcalenergy / hpluse) << " H/E " << (hcalenergy > ecalenergy)
2009  << " P/(H+E) " << (secpin / hpluse) << " HCAL ENE " << hcalenergy
2010  << " ECAL ENE " << ecalenergy << " secPIN " << secpin << " Algo Track "
2011  << trkRef->algo() << std::endl;
2012  remove_this_kf = true;
2013  }
2014  } else {
2015  if (isHoHE) {
2016  LOGDRESSED("PFEGammaAlgo") << "REJECTED TRACK FOR H/H+E, CLUSTER NOT IN SC"
2017  << " H/H+E " << (hcalenergy / hpluse) << " H/E " << (hcalenergy > ecalenergy)
2018  << " P/(H+E) " << (secpin / hpluse) << " HCAL ENE " << hcalenergy
2019  << " ECAL ENE " << ecalenergy << " secPIN " << secpin << " Algo Track "
2020  << trkRef->algo() << std::endl;
2021  remove_this_kf = true;
2022  }
2023  }
2024  }
2025  }
2026  }
2027  if (remove_this_kf) {
2028  secd_kf = RO.secondaryKFs.erase(secd_kf);
2029  if (secd_kf == RO.secondaryKFs.end())
2030  break;
2031  }
2032  }
2033 }
2034 
2035 bool PFEGammaAlgo::isPrimaryTrack(const reco::PFBlockElementTrack& KfEl, const reco::PFBlockElementGsfTrack& GsfEl) {
2036  bool isPrimary = false;
2037 
2038  const GsfPFRecTrackRef& gsfPfRef = GsfEl.GsftrackRefPF();
2039 
2040  if (gsfPfRef.isNonnull()) {
2041  const PFRecTrackRef& kfPfRef = KfEl.trackRefPF();
2042  PFRecTrackRef kfPfRef_fromGsf = (*gsfPfRef).kfPFRecTrackRef();
2043  if (kfPfRef.isNonnull() && kfPfRef_fromGsf.isNonnull()) {
2044  reco::TrackRef kfref = (*kfPfRef).trackRef();
2045  reco::TrackRef kfref_fromGsf = (*kfPfRef_fromGsf).trackRef();
2046  if (kfref.isNonnull() && kfref_fromGsf.isNonnull()) {
2047  if (kfref == kfref_fromGsf)
2048  isPrimary = true;
2049  }
2050  }
2051  }
2052 
2053  return isPrimary;
2054 }
Vector3DBase
Definition: Vector3DBase.h:8
muonTagProbeFilters_cff.matched
matched
Definition: muonTagProbeFilters_cff.py:62
reco::PFBlockElementTrack::trackType
bool trackType(TrackType trType) const override
Definition: PFBlockElementTrack.h:28
reco::PFCandidateEGammaExtra::setHadEnergy
void setHadEnergy(float val)
set the had energy. The cluster energies should be entered before
Definition: PFCandidateEGammaExtra.cc:74
PFEGammaAlgo::linkRefinableObjectPrimaryGSFTrackToECAL
void linkRefinableObjectPrimaryGSFTrackToECAL(ProtoEGObject &)
Definition: PFEGammaAlgo.cc:1222
reco::PFBlockElementGsfTrack::GsftrackRef
const reco::GsfTrackRef & GsftrackRef() const
Definition: PFBlockElementGsfTrack.h:47
reco::PFBlockElementGsfTrack::trackType
bool trackType(TrackType trType) const override
Definition: PFBlockElementGsfTrack.h:31
PFEGammaAlgo::mergeROsByAnyLink
void mergeROsByAnyLink(std::list< ProtoEGObject > &)
Definition: PFEGammaAlgo.cc:1079
PFEGammaAlgo::_splayedblock
std::vector< std::vector< FlaggedPtr< const reco::PFBlockElement > > > _splayedblock
Definition: PFEGammaAlgo.h:148
LOGERR
#define LOGERR(x)
Definition: PFEGammaAlgo.cc:45
mps_fire.i
i
Definition: mps_fire.py:355
PFEGammaAlgo::ProtoEGObject::ecal2ps
ClusterMap ecal2ps
Definition: PFEGammaAlgo.h:92
PFClusterWidthAlgo::pflowEtaWidth
double pflowEtaWidth() const
Definition: PFClusterWidthAlgo.h:15
PFEGammaAlgo::linkRefinableObjectConvSecondaryKFsToSecondaryKFs
void linkRefinableObjectConvSecondaryKFsToSecondaryKFs(ProtoEGObject &)
Definition: PFEGammaAlgo.cc:1397
PFEGammaAlgo::ProtoEGObject::parentSC
const PFSCElement * parentSC
Definition: PFEGammaAlgo.h:86
PFMuonAlgo.h
CommutativePairs::contains
bool contains(T const &a, T const &b) const
Definition: CommutativePairs.h:20
PFEGammaAlgo::ProtoEGObject::electronSeed
reco::ElectronSeedRef electronSeed
Definition: PFEGammaAlgo.h:87
f
double f[11][100]
Definition: MuScleFitUtils.cc:78
reco::PFCandidateEGammaExtra::setDeltaEta
void setDeltaEta(float val)
set the delta eta
Definition: PFCandidateEGammaExtra.cc:85
PFEGammaAlgo::EgammaObjects
Definition: PFEGammaAlgo.h:117
PFEGammaAlgo::cfg_
PFEGConfigInfo const & cfg_
Definition: PFEGammaAlgo.h:221
convertSQLitetoXML_cfg.output
output
Definition: convertSQLitetoXML_cfg.py:32
reco::SuperCluster
Definition: SuperCluster.h:18
PFEGammaAlgo::ProtoEGObject::brems
std::vector< PFBremElement const * > brems
Definition: PFEGammaAlgo.h:96
X
#define X(str)
Definition: MuonsGrabber.cc:38
edm::Ref::isAvailable
bool isAvailable() const
Definition: Ref.h:537
min
T min(T a, T b)
Definition: MathUtil.h:58
reco::Vertex::z
double z() const
z coordinate
Definition: Vertex.h:120
reco::PFCandidateEGammaExtra::addExtraNonConvTrack
void addExtraNonConvTrack(const reco::PFBlockRef &blk, const reco::PFBlockElementTrack &tkref)
track counting for electrons and photons
Definition: PFCandidateEGammaExtra.h:177
PFEGammaAlgo::channelStatus_
ESChannelStatus const & channelStatus_
Definition: PFEGammaAlgo.h:224
hcal
Definition: ConfigurationDatabase.cc:13
edm::Ref::isNull
bool isNull() const
Checks for null.
Definition: Ref.h:235
reco::PFBlockElementBrem::indTrajPoint
unsigned int indTrajPoint() const
Definition: PFBlockElementBrem.h:38
reco::CaloCluster::setCorrectedEnergy
void setCorrectedEnergy(double cenergy)
Definition: CaloCluster.h:137
reco::PFBlockElementSuperCluster
Cluster Element.
Definition: PFBlockElementSuperCluster.h:15
PFEGammaAlgo::linkRefinableObjectECALToSingleLegConv
void linkRefinableObjectECALToSingleLegConv(ProtoEGObject &)
Definition: PFEGammaAlgo.cc:1420
reco::PFBlock
Block of elements.
Definition: PFBlock.h:26
PFEGammaAlgo::ProtoEGObject::electronClusters
std::vector< const PFClusterElement * > electronClusters
Definition: PFEGammaAlgo.h:106
reco::SuperCluster::setEtaWidth
void setEtaWidth(double ew)
Definition: SuperCluster.h:74
EBDetId.h
EEDetId.h
reco::PFBlockElement::SC
Definition: PFBlockElement.h:41
reco::PFCandidateEGammaExtra::addConversionRef
void addConversionRef(const reco::ConversionRef &convref)
add Conversions from PF
Definition: PFCandidateEGammaExtra.cc:143
PFEGammaAlgo::PFEGConfigInfo::applyCrackCorrections
bool applyCrackCorrections
Definition: PFEGammaAlgo.h:112
reco::PFBlockElementTrack::trackRef
const reco::TrackRef & trackRef() const override
Definition: PFBlockElementTrack.h:49
edm::Ref::get
T const * get() const
Returns C++ pointer to the item.
Definition: Ref.h:232
HLT_2018_cff.distance
distance
Definition: HLT_2018_cff.py:6417
info
static const TGPicture * info(bool iBackgroundIsBlack)
Definition: FWCollectionSummaryWidget.cc:152
PFEGammaAlgo::ProtoEGObject::nBremsWithClusters
int nBremsWithClusters
Definition: PFEGammaAlgo.h:107
PFEGammaAlgo::gbrForests_
GBRForests const & gbrForests_
Definition: PFEGammaAlgo.h:135
reco::Vertex::position
const Point & position() const
position
Definition: Vertex.h:114
PFEGammaAlgo::ProtoEGObject::ecalclusters
std::vector< FlaggedPtr< const PFClusterElement > > ecalclusters
Definition: PFEGammaAlgo.h:91
spr::goodTrack
bool goodTrack(const reco::Track *pTrack, math::XYZPoint leadPV, trackSelectionParameters parameters, bool debug=false)
Definition: TrackSelection.cc:10
reco::PFBlockElement::HCAL
Definition: PFBlockElement.h:36
edm::Ptr::get
T const * get() const
Returns C++ pointer to the item.
Definition: Ptr.h:139
PFEGammaAlgo::ProtoEGObject::localMap
CommutativePairs< const reco::PFBlockElement * > localMap
Definition: PFEGammaAlgo.h:103
DDAxes::x
edm::RefVector< ConversionCollection >
ESDetId.h
reco::PFBlockElement::HGCAL
Definition: PFBlockElement.h:43
reco
fixed size matrix
Definition: AlignmentAlgorithmBase.h:45
PFEGammaAlgo::evaluateSingleLegMVA
float evaluateSingleLegMVA(const reco::PFBlockRef &blockref, const reco::Vertex &primaryVtx, unsigned int trackIndex)
Definition: PFEGammaAlgo.cc:459
hltPixelTracks_cff.chi2
chi2
Definition: hltPixelTracks_cff.py:25
reco::SuperCluster::setPreshowerEnergy
void setPreshowerEnergy(double preshowerEnergy)
Definition: SuperCluster.h:70
edm::Ptr::key
key_type key() const
Definition: Ptr.h:163
spr::find
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:19
edm::Ref::index
key_type index() const
Definition: Ref.h:253
groupFilesInBlocks.temp
list temp
Definition: groupFilesInBlocks.py:142
training_settings.idx
idx
Definition: training_settings.py:16
reco::PFCandidateEGammaExtra
Definition: PFCandidateEGammaExtra.h:25
reco::PFCandidateEGammaExtra::setSuperClusterRef
void setSuperClusterRef(reco::SuperClusterRef sc)
set reference to the corresponding supercluster
Definition: PFCandidateEGammaExtra.h:108
AlCaHLTBitMon_QueryRunRegistry.comp
comp
Definition: AlCaHLTBitMon_QueryRunRegistry.py:249
reco::CaloClusterPtr
edm::Ptr< CaloCluster > CaloClusterPtr
Definition: CaloClusterFwd.h:21
reco::PFBlockElementGsfTrack::Pout
const math::XYZTLorentzVector & Pout() const
Definition: PFBlockElementGsfTrack.h:55
PFRecHit.h
reco::PFCandidateEGammaExtra::setMVA
void setMVA(float val)
set the result (mostly for debugging)
Definition: PFCandidateEGammaExtra.cc:116
reco::PFBlockElementSuperCluster::Dump
void Dump(std::ostream &out=std::cout, const char *tab=" ") const override
print the object inside the element
Definition: PFBlockElementSuperCluster.cc:11
Mustache.h
edm::Ref< ElectronSeedCollection >
reco::PFBlockElementGsfTrack::positionAtECALEntrance
const math::XYZPointF & positionAtECALEntrance() const
Definition: PFBlockElementGsfTrack.h:50
PFLayer::ECAL_BARREL
Definition: PFLayer.h:33
reco::PFCandidateEGammaExtra::setEarlyBrem
void setEarlyBrem(float val)
set EarlyBrem
Definition: PFCandidateEGammaExtra.cc:69
reco::PFBlockElement::Type
Type
Definition: PFBlockElement.h:30
PFClusterWidthAlgo::pflowPhiWidth
double pflowPhiWidth() const
Definition: PFClusterWidthAlgo.h:14
PFEGammaAlgo::unlinkRefinableObjectKFandECALWithBadEoverP
void unlinkRefinableObjectKFandECALWithBadEoverP(ProtoEGObject &)
Definition: PFEGammaAlgo.cc:1890
PFLayer.h
reco::SuperCluster::addCluster
void addCluster(const CaloClusterPtr &r)
add reference to constituent BasicCluster
Definition: SuperCluster.h:122
edm::refToPtr
Ptr< typename C::value_type > refToPtr(Ref< C, typename C::value_type, refhelper::FindUsingAdvance< C, typename C::value_type > > const &ref)
Definition: RefToPtr.h:18
deltaR.h
reco::SuperCluster::setSeed
void setSeed(const CaloClusterPtr &r)
list of used xtals by DetId // now inherited by CaloCluster
Definition: SuperCluster.h:107
reco::SuperCluster::preshowerClustersBegin
CaloCluster_iterator preshowerClustersBegin() const
fist iterator over PreshowerCluster constituents
Definition: SuperCluster.h:92
patCandidatesForDimuonsSequences_cff.hcal
hcal
Definition: patCandidatesForDimuonsSequences_cff.py:37
GetRecoTauVFromDQM_MC_cff.cl
cl
Definition: GetRecoTauVFromDQM_MC_cff.py:38
BeamMonitor_cff.primaryVertex
primaryVertex
hltOfflineBeamSpot for HLTMON
Definition: BeamMonitor_cff.py:7
PFEGammaAlgo::ProtoEGObject::primaryGSFs
std::vector< PFGSFElement const * > primaryGSFs
Definition: PFEGammaAlgo.h:94
reco::PFCandidateEGammaExtra::setGsfElectronClusterRef
void setGsfElectronClusterRef(const reco::PFBlockRef &blk, const reco::PFBlockElementCluster &ref)
set gsf electron cluster ref
Definition: PFCandidateEGammaExtra.h:88
PFEnergyCalibration.h
Photon.h
reco::PFCandidateEGammaExtra::setKfTrackRef
void setKfTrackRef(const reco::TrackRef &ref)
set kf track reference
Definition: PFCandidateEGammaExtra.cc:51
ECAL
Definition: HCALResponse.h:21
alignCSCRings.s
s
Definition: alignCSCRings.py:92
test
Definition: SmallWORMDict.h:13
ctpps_dqm_sourceclient-live_cfg.test
test
Definition: ctpps_dqm_sourceclient-live_cfg.py:7
PFEGammaAlgo::primaryVertex_
reco::Vertex const & primaryVertex_
Definition: PFEGammaAlgo.h:222
PFEGammaAlgo::isPrimaryTrack
bool isPrimaryTrack(const reco::PFBlockElementTrack &KfEl, const reco::PFBlockElementGsfTrack &GsfEl)
Definition: PFEGammaAlgo.cc:2035
reco::PFCandidateEGammaExtra::setGsfTrackPout
void setGsfTrackPout(const math::XYZTLorentzVector &pout)
set the pout (not trivial to get from the GSF track)
Definition: PFCandidateEGammaExtra.cc:49
std::swap
void swap(edm::DataFrameContainer &lhs, edm::DataFrameContainer &rhs)
Definition: DataFrameContainer.h:209
reco::PFBlockElementTrack::trackRefPF
const PFRecTrackRef & trackRefPF() const override
Definition: PFBlockElementTrack.h:46
SiPixelRawToDigiRegional_cfi.deltaPhi
deltaPhi
Definition: SiPixelRawToDigiRegional_cfi.py:9
LOGDRESSED
#define LOGDRESSED(x)
Definition: PFEGammaAlgo.cc:46
vars
vars
Definition: DeepTauId.cc:163
PFEGammaAlgo::calculateEleMVA
float calculateEleMVA(const ProtoEGObject &, reco::PFCandidateEGammaExtra &) const
Definition: PFEGammaAlgo.cc:1601
reco::PFBlockElement::TRACK
Definition: PFBlockElement.h:32
PVValHelper::eta
Definition: PVValidationHelpers.h:69
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
PFEGammaAlgo::GBRForests::ele_
const std::unique_ptr< const GBRForest > ele_
Definition: PFEGammaAlgo.h:80
reco::Vertex::tracks_end
trackRef_iterator tracks_end() const
last iterator over tracks
Definition: Vertex.cc:73
PFEGammaAlgo::ClusterMap
std::unordered_map< PFClusterElement const *, std::vector< PFClusterElement const * > > ClusterMap
Definition: PFEGammaAlgo.h:72
CommutativePairs::concatenate
void concatenate(CommutativePairs< T > const &other)
Add the pairs from another CommutativePairs to this.
Definition: CommutativePairs.h:40
PFCluster.h
reco::CaloCluster::badHcalMarker
Definition: CaloCluster.h:50
PFEGammaAlgo::_currentblock
reco::PFBlockRef _currentblock
Definition: PFEGammaAlgo.h:144
reco::PFBlockElementBrem
Track Element.
Definition: PFBlockElementBrem.h:17
RecoTauValidation_cfi.posX
posX
Definition: RecoTauValidation_cfi.py:288
HCALHighEnergyHPDFilter_cfi.energy
energy
Definition: HCALHighEnergyHPDFilter_cfi.py:5
reco::PFCandidateEGammaExtra::setLateBrem
void setLateBrem(float val)
set LateBrem
Definition: PFCandidateEGammaExtra.cc:64
PFEGammaAlgo::fillExtraInfo
void fillExtraInfo(const ProtoEGObject &, reco::PFCandidateEGammaExtra &)
Definition: PFEGammaAlgo.cc:1765
VtxSmearedParameters_cfi.Phi
Phi
Definition: VtxSmearedParameters_cfi.py:112
PFClusterWidthAlgo::pflowSigmaEtaEta
double pflowSigmaEtaEta() const
Definition: PFClusterWidthAlgo.h:16
PFClusterWidthAlgo.h
HCAL
Definition: HCALResponse.h:21
b
double b
Definition: hdecay.h:118
PFEGammaAlgo::linkRefinableObjectBremTangentsToECAL
void linkRefinableObjectBremTangentsToECAL(ProtoEGObject &)
Definition: PFEGammaAlgo.cc:1332
reco::SuperCluster::setPreshowerEnergyPlane1
void setPreshowerEnergyPlane1(double preshowerEnergy1)
Definition: SuperCluster.h:71
PFEnergyCalibration::initAlphaGamma_ESplanes_fromDB
void initAlphaGamma_ESplanes_fromDB(const ESEEIntercalibConstants *esEEInterCalib)
Definition: PFEnergyCalibration.h:67
ClusterClusterMapping::overlap
static bool overlap(const reco::CaloCluster &sc1, const reco::CaloCluster &sc, float minfrac=0.01, bool debug=false)
Definition: ClusterClusterMapping.cc:4
PFLayer::Layer
Layer
layer definition
Definition: PFLayer.h:29
reco::PFCandidateEGammaExtra::setGsfTrackRef
void setGsfTrackRef(const reco::GsfTrackRef &ref)
set gsftrack reference
Definition: PFCandidateEGammaExtra.cc:47
MuonTCMETValueMapProducer_cff.nLayers
nLayers
Definition: MuonTCMETValueMapProducer_cff.py:38
reco::PFBlockElement::TrackType
TrackType
Definition: PFBlockElement.h:47
LOGWARN
#define LOGWARN(x)
Definition: PFEGammaAlgo.cc:44
PFEGammaAlgo::ProtoEGObject::hcalClusters
std::vector< PFClusterElement const * > hcalClusters
Definition: PFEGammaAlgo.h:102
PFEGammaAlgo::linkRefinableObjectGSFTracksToKFs
void linkRefinableObjectGSFTracksToKFs(ProtoEGObject &)
Definition: PFEGammaAlgo.cc:1159
reco::PFCandidateEGammaExtra::setSuperClusterPFECALRef
void setSuperClusterPFECALRef(reco::SuperClusterRef sc)
set reference to the corresponding supercluster
Definition: PFCandidateEGammaExtra.h:111
PFTrackAlgoTools.h
PFEGammaAlgo::ProtoEGObject::singleLegConversionMvaMap
KFValMap singleLegConversionMvaMap
Definition: PFEGammaAlgo.h:100
reco::PFBlockElementTrack::positionAtECALEntrance
const math::XYZPointF & positionAtECALEntrance() const
Definition: PFBlockElementTrack.h:42
math::XYZPoint
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
delPhi
static double delPhi(const double phi1, const double phi2)
Definition: FFTGenericScaleCalculator.cc:20
a
double a
Definition: hdecay.h:119
reco::SuperCluster::seed
const CaloClusterPtr & seed() const
seed BasicCluster
Definition: SuperCluster.h:77
reco::PFCandidateEGammaExtra::setSigmaEtaEta
void setSigmaEtaEta(float val)
set the sigmaetaeta
Definition: PFCandidateEGammaExtra.cc:80
PFEGammaAlgo::PFEGConfigInfo
Definition: PFEGammaAlgo.h:110
reco::PFBlockElement::GSF
Definition: PFBlockElement.h:37
SiStripPI::max
Definition: SiStripPayloadInspectorHelper.h:169
reco::Vertex::tracks_begin
trackRef_iterator tracks_begin() const
first iterator over tracks
Definition: Vertex.cc:71
PFEGammaAlgo::initializeProtoCands
void initializeProtoCands(std::list< ProtoEGObject > &)
Definition: PFEGammaAlgo.cc:655
reco::Vertex::x
double x() const
x coordinate
Definition: Vertex.h:116
edm::Ref::isNonnull
bool isNonnull() const
Checks for non-null.
Definition: Ref.h:238
PFTrackAlgoTools::isGoodForEGMPrimary
bool isGoodForEGMPrimary(const reco::TrackBase::TrackAlgorithm &)
Definition: PFTrackAlgoTools.cc:144
reco::CaloCluster::addHitAndFraction
void addHitAndFraction(DetId id, float fraction)
Definition: CaloCluster.h:203
PFEGammaAlgo::linkRefinableObjectPrimaryKFsToSecondaryKFs
void linkRefinableObjectPrimaryKFsToSecondaryKFs(ProtoEGObject &)
Definition: PFEGammaAlgo.cc:1191
StorageManager_cfg.e1
e1
Definition: StorageManager_cfg.py:16
math::XYZVector
XYZVectorD XYZVector
spatial vector with cartesian internal representation
Definition: Vector3D.h:31
FlaggedPtr
Definition: FlaggedPtr.h:5
reco::TrackExtraRef
edm::Ref< TrackExtraCollection > TrackExtraRef
persistent reference to a TrackExtra
Definition: TrackExtraFwd.h:16
reco::PFBlockElement::ECAL
Definition: PFBlockElement.h:35
cand
Definition: decayParser.h:34
RefToPtr.h
PFEGammaAlgo::GBRForests::singleLeg_
const std::unique_ptr< const GBRForest > singleLeg_
Definition: PFEGammaAlgo.h:81
p4
double p4[4]
Definition: TauolaWrapper.h:92
PFEGammaAlgo::removeOrLinkECALClustersToKFTracks
void removeOrLinkECALClustersToKFTracks()
Definition: PFEGammaAlgo.cc:1000
reco::PFBlockElement::T_FROM_GAMMACONV
Definition: PFBlockElement.h:47
edm::LogVerbatim
Definition: MessageLogger.h:297
reco::PFCandidateEGammaExtra::addSingleLegConvTrackRefMva
void addSingleLegConvTrackRefMva(const std::pair< reco::TrackRef, float > &trackrefmva)
add Single Leg Conversion TrackRef
Definition: PFCandidateEGammaExtra.cc:139
PFEGammaAlgo::attachPSClusters
int attachPSClusters(const PFClusterElement *, ClusterMap::mapped_type &)
Definition: PFEGammaAlgo.cc:901
ntupleEnum.gsf
gsf
Definition: ntupleEnum.py:48
groupFilesInBlocks.block
block
Definition: groupFilesInBlocks.py:150
PFEGammaAlgo::linkRefinableObjectKFTracksToECAL
void linkRefinableObjectKFTracksToECAL(ProtoEGObject &)
Definition: PFEGammaAlgo.cc:1286
CommutativePairs::insert
void insert(T const &a, T const &b)
Definition: CommutativePairs.h:17
HltBtagPostValidation_cff.c
c
Definition: HltBtagPostValidation_cff.py:31
PFEGammaAlgo::EEtoPSAssociation
reco::PFCluster::EEtoPSAssociation EEtoPSAssociation
Definition: PFEGammaAlgo.h:64
reco::SuperCluster::addPreshowerCluster
void addPreshowerCluster(const CaloClusterPtr &r)
add reference to constituent BasicCluster
Definition: SuperCluster.h:128
PFEGammaAlgo::linkRefinableObjectSecondaryKFsToECAL
void linkRefinableObjectSecondaryKFsToECAL(ProtoEGObject &)
Definition: PFEGammaAlgo.cc:1450
reco::PFBlock::LinkData
std::map< unsigned int, Link > LinkData
Definition: PFBlock.h:39
reco::PFBlock::LINKTEST_ALL
Definition: PFBlock.h:41
ESCondObjectContainer
Definition: ESCondObjectContainer.h:11
PFEGammaAlgo::GBRForests
Definition: PFEGammaAlgo.h:74
reco::PFBlockElement
Abstract base class for a PFBlock element (track, cluster...)
Definition: PFBlockElement.h:26
reco::PFBlockElementSuperCluster::superClusterRef
const SuperClusterRef & superClusterRef() const
Definition: PFBlockElementSuperCluster.h:36
cmsLHEtoEOSManager.l
l
Definition: cmsLHEtoEOSManager.py:193
reco::CaloCluster::position
const math::XYZPoint & position() const
cluster centroid position
Definition: CaloCluster.h:154
PFRecHitFraction.h
edm::Ptr< reco::PFCluster >
looper.cfg
cfg
Definition: looper.py:297
TtFullHadJetPartonMatch_cfi.maxDist
maxDist
Definition: TtFullHadJetPartonMatch_cfi.py:35
reco::SuperCluster::setPreshowerEnergyPlane2
void setPreshowerEnergyPlane2(double preshowerEnergy2)
Definition: SuperCluster.h:72
PFEGammaAlgo::ProtoEGObject::primaryKFs
std::vector< PFKFElement const * > primaryKFs
Definition: PFEGammaAlgo.h:95
reco::PFBlockElementGsfTrack
Track Element.
Definition: PFBlockElementGsfTrack.h:18
bookConverter.elements
elements
Definition: bookConverter.py:147
eostools.move
def move(src, dest)
Definition: eostools.py:511
std
Definition: JetResolutionObject.h:76
stdcomb::next_combination
bool next_combination(BidIt n_begin, BidIt n_end, BidIt r_begin, BidIt r_end)
Definition: combination.h:19
reco::Vertex::y
double y() const
y coordinate
Definition: Vertex.h:118
ClusterClusterMapping.h
PFMuonAlgo::isMuon
static bool isMuon(const reco::PFBlockElement &elt)
Definition: PFMuonAlgo.cc:56
PFEnergyResolution.h
PFEGammaAlgo::linkKFTrackToECAL
void linkKFTrackToECAL(PFKFElement const *, ProtoEGObject &)
Definition: PFEGammaAlgo.cc:1295
PFEGammaAlgo::ProtoEGObject::lateBrem
int lateBrem
Definition: PFEGammaAlgo.h:107
T
long double T
Definition: Basic3DVectorLD.h:48
math::XYZTLorentzVector
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:29
PFEGammaAlgo::linkRefinableObjectPrimaryGSFTrackToHCAL
void linkRefinableObjectPrimaryGSFTrackToHCAL(ProtoEGObject &)
Definition: PFEGammaAlgo.cc:1265
PFEGammaAlgo::ProtoEGObject::secondaryKFs
std::vector< PFKFElement const * > secondaryKFs
Definition: PFEGammaAlgo.h:99
PFTrackAlgoTools::isGoodForEGM
bool isGoodForEGM(const reco::TrackBase::TrackAlgorithm &)
Definition: PFTrackAlgoTools.cc:117
reco::SuperCluster::rawEnergy
double rawEnergy() const
raw uncorrected energy (sum of energies of component BasicClusters)
Definition: SuperCluster.h:58
relativeConstraints.empty
bool empty
Definition: relativeConstraints.py:46
reco::PFBlockElementCluster
Cluster Element.
Definition: PFBlockElementCluster.h:16
Exception
Definition: hltDiff.cc:246
PFEGammaAlgo::_currentlinks
reco::PFBlock::LinkData _currentlinks
Definition: PFEGammaAlgo.h:145
PFClusterWidthAlgo
Definition: PFClusterWidthAlgo.h:6
electronStore.links
links
Definition: electronStore.py:149
PFEGammaAlgo.h
reco::PFBlockElementTrack
Track Element.
Definition: PFBlockElementTrack.h:17
reco::SuperCluster::preshowerClustersEnd
CaloCluster_iterator preshowerClustersEnd() const
last iterator over PreshowerCluster constituents
Definition: SuperCluster.h:95
PFEGammaAlgo::ProtoEGObject::firstBrem
int firstBrem
Definition: PFEGammaAlgo.h:107
PFEGammaAlgo::isMuon
bool isMuon(const reco::PFBlockElement &)
Definition: PFEGammaAlgo.cc:505
reco::SuperCluster::setPhiWidth
void setPhiWidth(double pw)
Definition: SuperCluster.h:73
PFEGammaAlgo::PFEGConfigInfo::mvaConvCut
double mvaConvCut
Definition: PFEGammaAlgo.h:114
edm::Ptr::isNonnull
bool isNonnull() const
Checks for non-null.
Definition: Ptr.h:146
reco::PFBlockElementGsfTrack::GsftrackRefPF
const GsfPFRecTrackRef & GsftrackRefPF() const
Definition: PFBlockElementGsfTrack.h:44
reco::PFBlockElementCluster::clusterRef
const PFClusterRef & clusterRef() const override
Definition: PFBlockElementCluster.h:29
PFEGammaAlgo::ProtoEGObject::parentBlock
reco::PFBlockRef parentBlock
Definition: PFEGammaAlgo.h:85
dqm-mbProfile.log
log
Definition: dqm-mbProfile.py:17
PFEnergyCalibration::CalibratedEndcapPFClusterEnergies::clusterEnergy
double clusterEnergy
Definition: PFEnergyCalibration.h:50
reco::PFBlockElement::type
Type type() const
Definition: PFBlockElement.h:69
bsc_activity_cfg.ecal
ecal
Definition: bsc_activity_cfg.py:25
LOGVERB
#define LOGVERB(x)
Definition: PFEGammaAlgo.cc:43
electronIdCutBasedExt_cfi.sigmaEtaEta
sigmaEtaEta
Definition: electronIdCutBasedExt_cfi.py:123
PFEGammaAlgo::eetops_
reco::PFCluster::EEtoPSAssociation const & eetops_
Definition: PFEGammaAlgo.h:143
reco::PFCandidate
Particle reconstructed by the particle flow algorithm.
Definition: PFCandidate.h:40
mps_fire.result
result
Definition: mps_fire.py:303
reco::PFBlockElementGsfTrack::Dump
void Dump(std::ostream &out=std::cout, const char *tab=" ") const override
print the object inside the element
Definition: PFBlockElementGsfTrack.cc:30
ESEEIntercalibConstants
Definition: ESEEIntercalibConstants.h:7
reco::PFBlockElementGsfTrack::Pin
const math::XYZTLorentzVector & Pin() const
Definition: PFBlockElementGsfTrack.h:54
PFEGammaAlgo::unwrapSuperCluster
bool unwrapSuperCluster(const reco::PFBlockElementSuperCluster *, std::vector< FlaggedPtr< const PFClusterElement >> &, ClusterMap &)
Definition: PFEGammaAlgo.cc:776
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
reco::PFBlockElement::PS1
Definition: PFBlockElement.h:33
PFEGammaAlgo::operator()
EgammaObjects operator()(const reco::PFBlockRef &block)
Definition: PFEGammaAlgo.cc:527
reco::Vertex::index
unsigned int index
index type
Definition: Vertex.h:54
PFEnergyCalibration::calibrateEndcapClusterEnergies
CalibratedEndcapPFClusterEnergies calibrateEndcapClusterEnergies(reco::PFCluster const &eeCluster, std::vector< reco::PFCluster const * > const &psClusterPointers, ESChannelStatus const &channelStatus, bool applyCrackCorrections) const
Definition: PFEnergyCalibration.cc:154
reco::PFBlockElement::index
unsigned index() const
Definition: PFBlockElement.h:86
math::XYZPointF
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< float > > XYZPointF
point in space with cartesian internal representation
Definition: Point3D.h:10
PFLayer::ECAL_ENDCAP
Definition: PFLayer.h:32
docast
#define docast(x, y)
Definition: PFEGammaAlgo.cc:42
crabWrapper.key
key
Definition: crabWrapper.py:19
PFEGammaAlgo::PFEGammaAlgo
PFEGammaAlgo(const PFEGConfigInfo &, GBRForests const &gbrForests, EEtoPSAssociation const &eetops, ESEEIntercalibConstants const &esEEInterCalib, ESChannelStatus const &channelStatus, reco::Vertex const &primaryVertex)
Definition: PFEGammaAlgo.cc:445
PFEGammaAlgo::dumpCurrentRefinableObjects
void dumpCurrentRefinableObjects() const
Definition: PFEGammaAlgo.cc:931
PFEGammaAlgo::ProtoEGObject
Definition: PFEGammaAlgo.h:84
class-composition.parent
parent
Definition: class-composition.py:88
listHistos.trackPt
trackPt
Definition: listHistos.py:120
reco::CaloCluster::energy
double energy() const
cluster energy
Definition: CaloCluster.h:149
TauDecayModes.dec
dec
Definition: TauDecayModes.py:143
PV3DBase::phi
Geom::Phi< T > phi() const
Definition: PV3DBase.h:66
PFEGammaAlgo::buildRefinedSuperCluster
reco::SuperCluster buildRefinedSuperCluster(const ProtoEGObject &)
Definition: PFEGammaAlgo.cc:1786
reco::PFBlockElementSuperCluster::fromPFSuperCluster
bool fromPFSuperCluster() const
Definition: PFBlockElementSuperCluster.h:86
PFEGammaAlgo::fillPFCandidates
EgammaObjects fillPFCandidates(const std::list< ProtoEGObject > &)
Definition: PFEGammaAlgo.cc:1468
PFEGammaAlgo::unlinkRefinableObjectKFandECALMatchedToHCAL
void unlinkRefinableObjectKFandECALMatchedToHCAL(ProtoEGObject &, bool removeFreeECAL=false, bool removeSCECAL=false)
Definition: PFEGammaAlgo.cc:1963
reco::Vertex
Definition: Vertex.h:35
RecoTauValidation_cfi.posY
posY
Definition: RecoTauValidation_cfi.py:289
PFEGammaAlgo::PFEGConfigInfo::produceEGCandsWithNoSuperCluster
bool produceEGCandsWithNoSuperCluster
Definition: PFEGammaAlgo.h:113
deltaPhi.h
reco::PFBlockElementCluster::Dump
void Dump(std::ostream &out=std::cout, const char *tab=" ") const override
print the object inside the element
Definition: PFBlockElementCluster.cc:11
reco::PFBlockElement::PS2
Definition: PFBlockElement.h:34
edm::OwnVector< reco::PFBlockElement >
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
PFEGammaAlgo::thePFEnergyCalibration_
PFEnergyCalibration thePFEnergyCalibration_
Definition: PFEGammaAlgo.h:137
Generated for CMSSW Reference Manual by doxygen 1.8.16