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()->missingInnerHits();
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->missingInnerHits();
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.setVertex(RO.primaryKFs[0]->trackRef()->vertex());
1492  cand.addElementInBlock(_currentblock, RO.primaryKFs[0]->index());
1493  }
1494  if (!RO.primaryGSFs.empty()) {
1495  cand.setCharge(RO.primaryGSFs[0]->GsftrackRef()->chargeMode());
1496  xtra.setGsfTrackRef(RO.primaryGSFs[0]->GsftrackRef());
1497  cand.setGsfTrackRef(RO.primaryGSFs[0]->GsftrackRef());
1498  cand.setVertex(RO.primaryGSFs[0]->GsftrackRef()->vertex());
1499  cand.addElementInBlock(_currentblock, RO.primaryGSFs[0]->index());
1500  }
1501  if (RO.parentSC) {
1502  xtra.setSuperClusterPFECALRef(RO.parentSC->superClusterRef());
1503  // we'll set to the refined supercluster back up in the producer
1504  cand.setSuperClusterRef(RO.parentSC->superClusterRef());
1505  xtra.setSuperClusterRef(RO.parentSC->superClusterRef());
1506  cand.addElementInBlock(_currentblock, RO.parentSC->index());
1507  }
1508  // add brems
1509  for (const auto& brem : RO.brems) {
1510  cand.addElementInBlock(_currentblock, brem->index());
1511  }
1512  // add clusters and ps clusters
1513  for (const auto& ecal : RO.ecalclusters) {
1514  const PFClusterElement* clus = ecal.get();
1515  cand.addElementInBlock(_currentblock, clus->index());
1516  if (RO.ecal2ps.count(clus)) {
1517  for (auto& psclus : RO.ecal2ps.at(clus)) {
1518  cand.addElementInBlock(_currentblock, psclus->index());
1519  }
1520  }
1521  }
1522  // add secondary tracks
1523  for (const auto& kf : RO.secondaryKFs) {
1524  cand.addElementInBlock(_currentblock, kf->index());
1525  const reco::ConversionRefVector& convrefs = kf->convRefs();
1526  bool no_conv_ref = true;
1527  for (const auto& convref : convrefs) {
1528  if (convref.isNonnull() && convref.isAvailable()) {
1529  xtra.addConversionRef(convref);
1530  no_conv_ref = false;
1531  }
1532  }
1533  if (no_conv_ref) {
1534  //single leg conversions
1535 
1536  //look for stored mva value in map or else recompute
1537  const auto& mvavalmapped = RO.singleLegConversionMvaMap.find(kf);
1538  //FIXME: Abuse single mva value to store both provenance and single leg mva score
1539  //by storing 3.0 + mvaval
1540  float mvaval = (mvavalmapped != RO.singleLegConversionMvaMap.end()
1541  ? mvavalmapped->second
1542  : 3.0 + evaluateSingleLegMVA(_currentblock, primaryVertex_, kf->index()));
1543 
1544  xtra.addSingleLegConvTrackRefMva(std::make_pair(kf->trackRef(), mvaval));
1545  }
1546  }
1547 
1548  // build the refined supercluster from those clusters left in the cand
1549  output.refinedSuperClusters.push_back(buildRefinedSuperCluster(RO));
1550 
1551  //*TODO* cluster time is not reliable at the moment, so only use track timing
1552  float trkTime = 0, trkTimeErr = -1;
1553  if (!RO.primaryGSFs.empty() && RO.primaryGSFs[0]->isTimeValid()) {
1554  trkTime = RO.primaryGSFs[0]->time();
1555  trkTimeErr = RO.primaryGSFs[0]->timeError();
1556  } else if (!RO.primaryKFs.empty() && RO.primaryKFs[0]->isTimeValid()) {
1557  trkTime = RO.primaryKFs[0]->time();
1558  trkTimeErr = RO.primaryKFs[0]->timeError();
1559  }
1560  if (trkTimeErr >= 0) {
1561  cand.setTime(trkTime, trkTimeErr);
1562  }
1563 
1564  const reco::SuperCluster& the_sc = output.refinedSuperClusters.back();
1565  // with the refined SC in hand we build a naive candidate p4
1566  // and set the candidate ECAL position to either the barycenter of the
1567  // supercluster (if super-cluster present) or the seed of the
1568  // new SC generated by the EGAlgo
1569  const double scE = the_sc.energy();
1570  if (scE != 0.0) {
1571  const math::XYZPoint& seedPos = the_sc.seed()->position();
1572  math::XYZVector egDir = the_sc.position() - primaryVertex_.position();
1573  egDir = egDir.Unit();
1574  cand.setP4(math::XYZTLorentzVector(scE * egDir.x(), scE * egDir.y(), scE * egDir.z(), scE));
1575  math::XYZPointF ecalPOS_f(seedPos.x(), seedPos.y(), seedPos.z());
1576  cand.setPositionAtECALEntrance(ecalPOS_f);
1577  cand.setEcalEnergy(the_sc.rawEnergy(), the_sc.energy());
1578  } else if (cfg_.produceEGCandsWithNoSuperCluster && !RO.primaryGSFs.empty()) {
1579  const PFGSFElement* gsf = RO.primaryGSFs[0];
1580  const reco::GsfTrackRef& gref = gsf->GsftrackRef();
1581  math::XYZTLorentzVector p4(gref->pxMode(), gref->pyMode(), gref->pzMode(), gref->pMode());
1582  cand.setP4(p4);
1583  cand.setPositionAtECALEntrance(gsf->positionAtECALEntrance());
1584  } else if (cfg_.produceEGCandsWithNoSuperCluster && !RO.primaryKFs.empty()) {
1585  const PFKFElement* kf = RO.primaryKFs[0];
1586  reco::TrackRef kref = RO.primaryKFs[0]->trackRef();
1587  math::XYZTLorentzVector p4(kref->px(), kref->py(), kref->pz(), kref->p());
1588  cand.setP4(p4);
1589  cand.setPositionAtECALEntrance(kf->positionAtECALEntrance());
1590  }
1591  const float eleMVAValue = calculateEleMVA(RO, xtra);
1592  fillExtraInfo(RO, xtra);
1593  //std::cout << "PFEG eleMVA: " << eleMVAValue << std::endl;
1594  xtra.setMVA(eleMVAValue);
1595  cand.set_mva_e_pi(eleMVAValue);
1596  output.candidates.push_back(cand);
1597  output.candidateExtras.push_back(xtra);
1598  }
1599 
1600  return output;
1601 }
1602 
1603 float PFEGammaAlgo::calculateEleMVA(const PFEGammaAlgo::ProtoEGObject& ro, reco::PFCandidateEGammaExtra& xtra) const {
1604  if (ro.primaryGSFs.empty()) {
1605  return -2.0f;
1606  }
1607  const PFGSFElement* gsfElement = ro.primaryGSFs.front();
1608  const PFKFElement* kfElement = nullptr;
1609  if (!ro.primaryKFs.empty()) {
1610  kfElement = ro.primaryKFs.front();
1611  }
1612  auto const& refGsf = gsfElement->GsftrackRef();
1613  reco::TrackRef refKf;
1614  constexpr float mEl = 0.000511;
1615  const double eInGsf = std::hypot(refGsf->pMode(), mEl);
1616  double dEtGsfEcal = 1e6;
1617  double sigmaEtaEta = 1e-14;
1618  const double eneHcalGsf =
1619  std::accumulate(ro.hcalClusters.begin(), ro.hcalClusters.end(), 0.0, [](const double a, auto const& b) {
1620  return a + b->clusterRef()->energy();
1621  });
1622  if (!ro.primaryKFs.empty()) {
1623  refKf = ro.primaryKFs.front()->trackRef();
1624  }
1625  const double eOutGsf = gsfElement->Pout().t();
1626  const double etaOutGsf = gsfElement->positionAtECALEntrance().eta();
1627  double firstEcalGsfEnergy{0.0};
1628  double otherEcalGsfEnergy{0.0};
1629  double ecalBremEnergy{0.0};
1630  //shower shape of cluster closest to gsf track
1631  std::vector<const reco::PFCluster*> gsfCluster;
1632  for (const auto& ecal : ro.ecalclusters) {
1633  const double cenergy = ecal->clusterRef()->correctedEnergy();
1634  bool hasgsf = ro.localMap.contains(gsfElement, ecal.get());
1635  bool haskf = ro.localMap.contains(kfElement, ecal.get());
1636  bool hasbrem = false;
1637  for (const auto& brem : ro.brems) {
1638  if (ro.localMap.contains(brem, ecal.get())) {
1639  hasbrem = true;
1640  }
1641  }
1642  if (hasbrem && ecal.get() != ro.electronClusters[0]) {
1643  ecalBremEnergy += cenergy;
1644  }
1645  if (!hasbrem && ecal.get() != ro.electronClusters[0]) {
1646  if (hasgsf)
1647  otherEcalGsfEnergy += cenergy;
1648  if (haskf)
1649  ecalBremEnergy += cenergy; // from conv. brem!
1650  if (!(hasgsf || haskf))
1651  otherEcalGsfEnergy += cenergy; // stuff from SC
1652  }
1653  }
1654 
1655  if (ro.electronClusters[0]) {
1656  reco::PFClusterRef cref = ro.electronClusters[0]->clusterRef();
1657  xtra.setGsfElectronClusterRef(_currentblock, *(ro.electronClusters[0]));
1658  firstEcalGsfEnergy = cref->correctedEnergy();
1659  dEtGsfEcal = cref->positionREP().eta() - etaOutGsf;
1660  gsfCluster.push_back(cref.get());
1661  PFClusterWidthAlgo pfwidth(gsfCluster);
1662  sigmaEtaEta = pfwidth.pflowSigmaEtaEta();
1663  }
1664 
1665  // brem sequence information
1666  float firstBrem{-1.0f};
1667  float earlyBrem{-1.0f};
1668  float lateBrem{-1.0f};
1669  if (ro.nBremsWithClusters > 0) {
1670  firstBrem = ro.firstBrem;
1671  earlyBrem = ro.firstBrem < 4 ? 1.0f : 0.0f;
1672  lateBrem = ro.lateBrem == 1 ? 1.0f : 0.0f;
1673  }
1674  xtra.setEarlyBrem(earlyBrem);
1675  xtra.setLateBrem(lateBrem);
1676  if (firstEcalGsfEnergy > 0.0) {
1677  if (refGsf.isNonnull()) {
1678  xtra.setGsfTrackPout(gsfElement->Pout());
1679  // normalization observables
1680  const float ptGsf = refGsf->ptMode();
1681  const float etaGsf = refGsf->etaMode();
1682  // tracking observables
1683  const double ptModeErrorGsf = refGsf->ptModeError();
1684  float ptModeErrOverPtGsf = (ptModeErrorGsf > 0. ? ptModeErrorGsf / ptGsf : 1.0);
1685  float chi2Gsf = refGsf->normalizedChi2();
1686  float dPtOverPtGsf = (ptGsf - gsfElement->Pout().pt()) / ptGsf;
1687  // kalman filter vars
1688  float nHitKf = refKf.isNonnull() ? refKf->hitPattern().trackerLayersWithMeasurement() : 0;
1689  float chi2Kf = refKf.isNonnull() ? refKf->normalizedChi2() : -0.01;
1690 
1691  //tracker + calorimetry observables
1692  float eTotPinMode = (firstEcalGsfEnergy + otherEcalGsfEnergy + ecalBremEnergy) / eInGsf;
1693  float eGsfPoutMode = firstEcalGsfEnergy / eOutGsf;
1694  float eTotBremPinPoutMode = (ecalBremEnergy + otherEcalGsfEnergy) / (eInGsf - eOutGsf);
1695  float dEtaGsfEcalClust = std::abs(dEtGsfEcal);
1696  float logSigmaEtaEta = std::log(sigmaEtaEta);
1697  float hOverHe = eneHcalGsf / (eneHcalGsf + firstEcalGsfEnergy);
1698 
1699  xtra.setDeltaEta(dEtaGsfEcalClust);
1700  xtra.setSigmaEtaEta(sigmaEtaEta);
1701  xtra.setHadEnergy(eneHcalGsf);
1702 
1703  // Apply bounds to variables and calculate MVA
1704  dPtOverPtGsf = std::clamp(dPtOverPtGsf, -0.2f, 1.0f);
1705  ptModeErrOverPtGsf = std::min(ptModeErrOverPtGsf, 0.3f);
1706  chi2Gsf = std::min(chi2Gsf, 10.0f);
1707  chi2Kf = std::min(chi2Kf, 10.0f);
1708  eTotPinMode = std::clamp(eTotPinMode, 0.0f, 5.0f);
1709  eGsfPoutMode = std::clamp(eGsfPoutMode, 0.0f, 5.0f);
1710  eTotBremPinPoutMode = std::clamp(eTotBremPinPoutMode, 0.0f, 5.0f);
1711  dEtaGsfEcalClust = std::min(dEtaGsfEcalClust, 0.1f);
1712  logSigmaEtaEta = std::max(logSigmaEtaEta, -14.0f);
1713 
1714  // not used for moment, weird behavior of variable
1715  //float dPtOverPtKf = refKf.isNonnull() ? (refKf->pt() - refKf->outerPt())/refKf->pt() : -0.01;
1716  //dPtOverPtKf = std::clamp(dPtOverPtKf,-0.2f, 1.0f);
1717 
1718  /*
1719  * To be used for debugging:
1720  * pretty-print the PFEgamma electron MVA input variables
1721  *
1722  * std::cout << " **** PFEG BDT observables ****" << endl;
1723  * std::cout << " < Normalization > " << endl;
1724  * std::cout << " ptGsf " << ptGsf << " Pin " << eInGsf
1725  * << " Pout " << eOutGsf << " etaGsf " << etaGsf << endl;
1726  * std::cout << " < PureTracking > " << endl;
1727  * std::cout << " ptModeErrOverPtGsf " << ptModeErrOverPtGsf
1728  * << " dPtOverPtGsf " << dPtOverPtGsf
1729  * << " chi2Gsf " << chi2Gsf
1730  * << " nhit_gsf " << nhit_gsf
1731  * << " dPtOverPtKf " << dPtOverPtKf
1732  * << " chi2Kf " << chi2Kf
1733  * << " nHitKf " << nHitKf << endl;
1734  * std::cout << " < track-ecal-hcal-ps " << endl;
1735  * std::cout << " eTotPinMode " << eTotPinMode
1736  * << " eGsfPoutMode " << eGsfPoutMode
1737  * << " eTotBremPinPoutMode " << eTotBremPinPoutMode
1738  * << " dEtaGsfEcalClust " << dEtaGsfEcalClust
1739  * << " logSigmaEtaEta " << logSigmaEtaEta
1740  * << " hOverHe " << hOverHe << " Hcal energy " << eneHcalGsf
1741  * << " lateBrem " << lateBrem
1742  * << " firstBrem " << firstBrem << endl;
1743  */
1744 
1745  float vars[] = {std::log(ptGsf),
1746  etaGsf,
1747  ptModeErrOverPtGsf,
1748  dPtOverPtGsf,
1749  chi2Gsf,
1750  nHitKf,
1751  chi2Kf,
1752  eTotPinMode,
1753  eGsfPoutMode,
1754  eTotBremPinPoutMode,
1755  dEtaGsfEcalClust,
1756  logSigmaEtaEta,
1757  hOverHe,
1758  lateBrem,
1759  firstBrem};
1760 
1761  return gbrForests_.ele_->GetAdaBoostClassifier(vars);
1762  }
1763  }
1764  return -2.0f;
1765 }
1766 
1767 void PFEGammaAlgo::fillExtraInfo(const ProtoEGObject& RO, reco::PFCandidateEGammaExtra& xtra) {
1768  // add tracks associated to clusters that are not T_FROM_GAMMACONV
1769  // info about single-leg convs is already save, so just veto in loops
1770  auto KFbegin = _splayedblock[reco::PFBlockElement::TRACK].begin();
1771  auto KFend = _splayedblock[reco::PFBlockElement::TRACK].end();
1772  for (auto& ecal : RO.ecalclusters) {
1773  NotCloserToOther<reco::PFBlockElement::ECAL, reco::PFBlockElement::TRACK, true> ECALToTracks(_currentblock,
1774  ecal.get());
1775  auto notmatchedkf = std::partition(KFbegin, KFend, ECALToTracks);
1776  auto notconvkf = std::partition(KFbegin, notmatchedkf, [](auto const& x) { return x->trackType(ConvType); });
1777  // go through non-conv-identified kfs and check MVA to add conversions
1778  for (auto kf = notconvkf; kf != notmatchedkf; ++kf) {
1779  const reco::PFBlockElementTrack* elemaskf = docast(const reco::PFBlockElementTrack*, kf->get());
1780  xtra.addExtraNonConvTrack(_currentblock, *elemaskf);
1781  }
1782  }
1783 }
1784 
1785 // currently stolen from PFECALSuperClusterAlgo, we should
1786 // try to factor this correctly since the operation is the same in
1787 // both places...
1788 reco::SuperCluster PFEGammaAlgo::buildRefinedSuperCluster(const PFEGammaAlgo::ProtoEGObject& RO) {
1789  if (RO.ecalclusters.empty()) {
1790  return reco::SuperCluster(0.0, math::XYZPoint(0, 0, 0));
1791  }
1792 
1793  bool isEE = false;
1794  // need the vector of raw pointers for a PF width class
1795  std::vector<const reco::PFCluster*> bare_ptrs;
1796  // calculate necessary parameters and build the SC
1797  double posX(0), posY(0), posZ(0), rawSCEnergy(0), corrSCEnergy(0), corrPSEnergy(0), ps1_energy(0.0), ps2_energy(0.0);
1798  for (auto& clus : RO.ecalclusters) {
1799  double ePS1 = 0;
1800  double ePS2 = 0;
1801  isEE = PFLayer::ECAL_ENDCAP == clus->clusterRef()->layer();
1802  auto clusptr = edm::refToPtr<reco::PFClusterCollection>(clus->clusterRef());
1803  bare_ptrs.push_back(clusptr.get());
1804 
1805  const double cluseraw = clusptr->energy();
1806  double cluscalibe = clusptr->correctedEnergy();
1807  const math::XYZPoint& cluspos = clusptr->position();
1808  posX += cluseraw * cluspos.X();
1809  posY += cluseraw * cluspos.Y();
1810  posZ += cluseraw * cluspos.Z();
1811  // update EE calibrated super cluster energies
1812  if (isEE && RO.ecal2ps.count(clus.get())) {
1813  const auto& psclusters = RO.ecal2ps.at(clus.get());
1814 
1815  std::vector<reco::PFCluster const*> psClusterPointers;
1816  psClusterPointers.reserve(psclusters.size());
1817  for (auto const& psc : psclusters) {
1818  psClusterPointers.push_back(psc->clusterRef().get());
1819  }
1820  auto calibratedEnergies = thePFEnergyCalibration_.calibrateEndcapClusterEnergies(
1821  *clusptr, psClusterPointers, channelStatus_, cfg_.applyCrackCorrections);
1822  cluscalibe = calibratedEnergies.clusterEnergy;
1823  ePS1 = calibratedEnergies.ps1Energy;
1824  ePS2 = calibratedEnergies.ps2Energy;
1825  }
1826  if (ePS1 == -1.)
1827  ePS1 = 0;
1828  if (ePS2 == -1.)
1829  ePS2 = 0;
1830 
1831  rawSCEnergy += cluseraw;
1832  corrSCEnergy += cluscalibe;
1833  ps1_energy += ePS1;
1834  ps2_energy += ePS2;
1835  corrPSEnergy += ePS1 + ePS2;
1836  }
1837  posX /= rawSCEnergy;
1838  posY /= rawSCEnergy;
1839  posZ /= rawSCEnergy;
1840 
1841  // now build the supercluster
1842  reco::SuperCluster new_sc(corrSCEnergy, math::XYZPoint(posX, posY, posZ));
1843 
1844  auto clusptr = edm::refToPtr<reco::PFClusterCollection>(RO.ecalclusters.front()->clusterRef());
1845  new_sc.setCorrectedEnergy(corrSCEnergy);
1846  new_sc.setSeed(clusptr);
1847  new_sc.setPreshowerEnergyPlane1(ps1_energy);
1848  new_sc.setPreshowerEnergyPlane2(ps2_energy);
1849  new_sc.setPreshowerEnergy(corrPSEnergy);
1850  for (const auto& clus : RO.ecalclusters) {
1851  clusptr = edm::refToPtr<reco::PFClusterCollection>(clus->clusterRef());
1852  new_sc.addCluster(clusptr);
1853  auto& hits_and_fractions = clusptr->hitsAndFractions();
1854  for (auto& hit_and_fraction : hits_and_fractions) {
1855  new_sc.addHitAndFraction(hit_and_fraction.first, hit_and_fraction.second);
1856  }
1857  // put the preshower stuff back in later
1858  if (RO.ecal2ps.count(clus.get())) {
1859  const auto& cluspsassociation = RO.ecal2ps.at(clus.get());
1860  // EE rechits should be uniquely matched to sets of pre-shower
1861  // clusters at this point, so we throw an exception if otherwise
1862  // now wrapped in EDM debug flags
1863  for (const auto& pscluselem : cluspsassociation) {
1864  edm::Ptr<reco::PFCluster> psclus = edm::refToPtr<reco::PFClusterCollection>(pscluselem->clusterRef());
1865 #ifdef PFFLOW_DEBUG
1866  auto found_pscluster =
1867  std::find(new_sc.preshowerClustersBegin(), new_sc.preshowerClustersEnd(), reco::CaloClusterPtr(psclus));
1868  if (found_pscluster == new_sc.preshowerClustersEnd()) {
1869 #endif
1870  new_sc.addPreshowerCluster(psclus);
1871 #ifdef PFFLOW_DEBUG
1872  } else {
1873  throw cms::Exception("PFECALSuperClusterAlgo::buildSuperCluster")
1874  << "Found a PS cluster matched to more than one EE cluster!" << std::endl
1875  << std::hex << psclus.get() << " == " << found_pscluster->get() << std::dec << std::endl;
1876  }
1877 #endif
1878  }
1879  }
1880  }
1881 
1882  // calculate linearly weighted cluster widths
1883  PFClusterWidthAlgo pfwidth(bare_ptrs);
1884  new_sc.setEtaWidth(pfwidth.pflowEtaWidth());
1885  new_sc.setPhiWidth(pfwidth.pflowPhiWidth());
1886 
1887  // cache the value of the raw energy
1888  new_sc.rawEnergy();
1889 
1890  return new_sc;
1891 }
1892 
1893 void PFEGammaAlgo::unlinkRefinableObjectKFandECALWithBadEoverP(ProtoEGObject& RO) {
1894  // this only means something for ROs with a primary GSF track
1895  if (RO.primaryGSFs.empty())
1896  return;
1897  // need energy sums to tell if we've added crap or not
1898  const double Pin_gsf = RO.primaryGSFs.front()->GsftrackRef()->pMode();
1899  const double gsfOuterEta = RO.primaryGSFs.front()->positionAtECALEntrance().Eta();
1900  double tot_ecal = 0.0;
1901  std::vector<double> min_brem_dists;
1902  std::vector<double> closest_brem_eta;
1903  // first get the total ecal energy (we should replace this with a cache)
1904  for (const auto& ecal : RO.ecalclusters) {
1905  tot_ecal += ecal->clusterRef()->correctedEnergy();
1906  // we also need to look at the minimum distance to brems
1907  // since energetic brems will be closer to the brem than the track
1908  double min_brem_dist = 5000.0;
1909  double eta = -999.0;
1910  for (const auto& brem : RO.brems) {
1911  const float dist = _currentblock->dist(brem->index(), ecal->index(), _currentlinks, reco::PFBlock::LINKTEST_ALL);
1912  if (dist < min_brem_dist && dist != -1.0f) {
1913  min_brem_dist = dist;
1914  eta = brem->positionAtECALEntrance().Eta();
1915  }
1916  }
1917  min_brem_dists.push_back(min_brem_dist);
1918  closest_brem_eta.push_back(eta);
1919  }
1920 
1921  // loop through the ECAL clusters and remove ECAL clusters matched to
1922  // secondary track either in *or* out of the SC if the E/pin is bad
1923  for (auto secd_kf = RO.secondaryKFs.begin(); secd_kf != RO.secondaryKFs.end(); ++secd_kf) {
1924  reco::TrackRef trkRef = (*secd_kf)->trackRef();
1925  const float secpin = (*secd_kf)->trackRef()->p();
1926  bool remove_this_kf = false;
1927  for (auto ecal = RO.ecalclusters.begin(); ecal != RO.ecalclusters.end(); ++ecal) {
1928  size_t bremidx = std::distance(RO.ecalclusters.begin(), ecal);
1929  const float minbremdist = min_brem_dists[bremidx];
1930  const double ecalenergy = (*ecal)->clusterRef()->correctedEnergy();
1931  const double Epin = ecalenergy / secpin;
1932  const double detaGsf = std::abs(gsfOuterEta - (*ecal)->clusterRef()->positionREP().Eta());
1933  const double detaBrem = std::abs(closest_brem_eta[bremidx] - (*ecal)->clusterRef()->positionREP().Eta());
1934 
1935  bool kf_matched = RO.localMap.contains(ecal->get(), *secd_kf);
1936 
1937  const float tkdist =
1938  _currentblock->dist((*secd_kf)->index(), (*ecal)->index(), _currentlinks, reco::PFBlock::LINKTEST_ALL);
1939 
1940  // do not reject this track if it is closer to a brem than the
1941  // secondary track, or if it lies in the delta-eta plane with the
1942  // gsf track or if it is in the dEta plane with the brems
1943  if (Epin > 3 && kf_matched && tkdist != -1.0f && tkdist < minbremdist && detaGsf > 0.05 && detaBrem > 0.015) {
1944  double res_with = std::abs((tot_ecal - Pin_gsf) / Pin_gsf);
1945  double res_without = std::abs((tot_ecal - ecalenergy - Pin_gsf) / Pin_gsf);
1946  if (res_without < res_with) {
1947  LOGDRESSED("PFEGammaAlgo") << " REJECTED_RES totenergy " << tot_ecal << " Pin_gsf " << Pin_gsf
1948  << " cluster to secondary " << ecalenergy << " res_with " << res_with
1949  << " res_without " << res_without << std::endl;
1950  tot_ecal -= ecalenergy;
1951  remove_this_kf = true;
1952  ecal = RO.ecalclusters.erase(ecal);
1953  if (ecal == RO.ecalclusters.end())
1954  break;
1955  }
1956  }
1957  }
1958  if (remove_this_kf) {
1959  secd_kf = RO.secondaryKFs.erase(secd_kf);
1960  if (secd_kf == RO.secondaryKFs.end())
1961  break;
1962  }
1963  }
1964 }
1965 
1966 void PFEGammaAlgo::unlinkRefinableObjectKFandECALMatchedToHCAL(ProtoEGObject& RO,
1967  bool removeFreeECAL,
1968  bool removeSCEcal) {
1969  std::vector<bool> cluster_in_sc;
1970  auto ecal_begin = RO.ecalclusters.begin();
1971  auto ecal_end = RO.ecalclusters.end();
1972  auto hcal_begin = _splayedblock[reco::PFBlockElement::HCAL].begin();
1973  auto hcal_end = _splayedblock[reco::PFBlockElement::HCAL].end();
1974  for (auto secd_kf = RO.secondaryKFs.begin(); secd_kf != RO.secondaryKFs.end(); ++secd_kf) {
1975  bool remove_this_kf = false;
1976  NotCloserToOther<reco::PFBlockElement::TRACK, reco::PFBlockElement::HCAL> tracksToHCALs(_currentblock, *secd_kf);
1977  reco::TrackRef trkRef = (*secd_kf)->trackRef();
1978 
1979  bool goodTrack = PFTrackAlgoTools::isGoodForEGM(trkRef->algo());
1980  const float secpin = trkRef->p();
1981 
1982  for (auto ecal = ecal_begin; ecal != ecal_end; ++ecal) {
1983  const double ecalenergy = (*ecal)->clusterRef()->correctedEnergy();
1984  // first check if the cluster is in the SC (use dist calc for fastness)
1985  const size_t clus_idx = std::distance(ecal_begin, ecal);
1986  if (cluster_in_sc.size() < clus_idx + 1) {
1987  float dist = -1.0f;
1988  if (RO.parentSC) {
1989  dist = _currentblock->dist((*secd_kf)->index(), (*ecal)->index(), _currentlinks, reco::PFBlock::LINKTEST_ALL);
1990  }
1991  cluster_in_sc.push_back(dist != -1.0f);
1992  }
1993 
1994  // if we've found a secondary KF that matches this ecal cluster
1995  // now we see if it is matched to HCAL
1996  // if it is matched to an HCAL cluster we take different
1997  // actions if the cluster was in an SC or not
1998  if (RO.localMap.contains(ecal->get(), *secd_kf)) {
1999  auto hcal_matched = std::partition(hcal_begin, hcal_end, tracksToHCALs);
2000  for (auto hcalclus = hcal_begin; hcalclus != hcal_matched; ++hcalclus) {
2001  const reco::PFBlockElementCluster* clusthcal =
2002  dynamic_cast<const reco::PFBlockElementCluster*>(hcalclus->get());
2003  const double hcalenergy = clusthcal->clusterRef()->energy();
2004  const double hpluse = ecalenergy + hcalenergy;
2005  const bool isHoHE = ((hcalenergy / hpluse) > 0.1 && goodTrack);
2006  const bool isHoE = (hcalenergy > ecalenergy);
2007  const bool isPoHE = (secpin > hpluse);
2008  if (cluster_in_sc[clus_idx]) {
2009  if (isHoE || isPoHE) {
2010  LOGDRESSED("PFEGammaAlgo") << "REJECTED TRACK FOR H/E or P/(H+E), CLUSTER IN SC"
2011  << " H/H+E " << (hcalenergy / hpluse) << " H/E " << (hcalenergy > ecalenergy)
2012  << " P/(H+E) " << (secpin / hpluse) << " HCAL ENE " << hcalenergy
2013  << " ECAL ENE " << ecalenergy << " secPIN " << secpin << " Algo Track "
2014  << trkRef->algo() << std::endl;
2015  remove_this_kf = true;
2016  }
2017  } else {
2018  if (isHoHE) {
2019  LOGDRESSED("PFEGammaAlgo") << "REJECTED TRACK FOR H/H+E, CLUSTER NOT IN SC"
2020  << " H/H+E " << (hcalenergy / hpluse) << " H/E " << (hcalenergy > ecalenergy)
2021  << " P/(H+E) " << (secpin / hpluse) << " HCAL ENE " << hcalenergy
2022  << " ECAL ENE " << ecalenergy << " secPIN " << secpin << " Algo Track "
2023  << trkRef->algo() << std::endl;
2024  remove_this_kf = true;
2025  }
2026  }
2027  }
2028  }
2029  }
2030  if (remove_this_kf) {
2031  secd_kf = RO.secondaryKFs.erase(secd_kf);
2032  if (secd_kf == RO.secondaryKFs.end())
2033  break;
2034  }
2035  }
2036 }
2037 
2038 bool PFEGammaAlgo::isPrimaryTrack(const reco::PFBlockElementTrack& KfEl, const reco::PFBlockElementGsfTrack& GsfEl) {
2039  bool isPrimary = false;
2040 
2041  const GsfPFRecTrackRef& gsfPfRef = GsfEl.GsftrackRefPF();
2042 
2043  if (gsfPfRef.isNonnull()) {
2044  const PFRecTrackRef& kfPfRef = KfEl.trackRefPF();
2045  PFRecTrackRef kfPfRef_fromGsf = (*gsfPfRef).kfPFRecTrackRef();
2046  if (kfPfRef.isNonnull() && kfPfRef_fromGsf.isNonnull()) {
2047  reco::TrackRef kfref = (*kfPfRef).trackRef();
2048  reco::TrackRef kfref_fromGsf = (*kfPfRef_fromGsf).trackRef();
2049  if (kfref.isNonnull() && kfref_fromGsf.isNonnull()) {
2050  if (kfref == kfref_fromGsf)
2051  isPrimary = true;
2052  }
2053  }
2054  }
2055 
2056  return isPrimary;
2057 }
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:428
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:72
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
deltaPhi.h
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
info
static const TGPicture * info(bool iBackgroundIsBlack)
Definition: FWCollectionSummaryWidget.cc:153
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
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
heavyIonCSV_trainingSettings.idx
idx
Definition: heavyIonCSV_trainingSettings.py:5
PFEGammaAlgo::unlinkRefinableObjectKFandECALWithBadEoverP
void unlinkRefinableObjectKFandECALWithBadEoverP(ProtoEGObject &)
Definition: PFEGammaAlgo.cc:1893
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
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:2038
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:159
PFEGammaAlgo::calculateEleMVA
float calculateEleMVA(const ProtoEGObject &, reco::PFCandidateEGammaExtra &) const
Definition: PFEGammaAlgo.cc:1603
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:1767
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
deltaR.h
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
phase2tkutil::isPrimary
bool isPrimary(const SimTrack &simTrk, const PSimHit *simHit)
Definition: TrackerPhase2ValidationUtil.cc:2
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:32
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
trackerHitRTTI::vector
Definition: trackerHitRTTI.h:21
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:204
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:48
PFEnergyResolution.h
PFEGammaAlgo::linkKFTrackToECAL
void linkKFTrackToECAL(PFKFElement const *, ProtoEGObject &)
Definition: PFEGammaAlgo.cc:1295
PFEGammaAlgo::ProtoEGObject::lateBrem
int lateBrem
Definition: PFEGammaAlgo.h:107
edm::LogVerbatim
Log< level::Info, true > LogVerbatim
Definition: MessageLogger.h:128
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:41
mps_fire.result
result
Definition: mps_fire.py:311
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
HLT_FULL_cff.distance
distance
Definition: HLT_FULL_cff.py:7799
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:1788
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:1966
reco::Vertex
Definition: Vertex.h:35
RecoTauValidation_cfi.posY
posY
Definition: RecoTauValidation_cfi.py:289
PFEGammaAlgo::PFEGConfigInfo::produceEGCandsWithNoSuperCluster
bool produceEGCandsWithNoSuperCluster
Definition: PFEGammaAlgo.h:113
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