CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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/PFRecHitFwd.h"
5 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
6 #include "DataFormats/ParticleFlowReco/interface/PFClusterFwd.h"
7 #include "RecoParticleFlow/PFClusterTools/interface/ClusterClusterMapping.h"
8 #include "DataFormats/ParticleFlowReco/interface/PFRecHit.h"
9 #include "DataFormats/EcalDetId/interface/EBDetId.h"
10 #include "DataFormats/EcalDetId/interface/EEDetId.h"
11 #include "DataFormats/EgammaCandidates/interface/Photon.h"
12 #include "RecoParticleFlow/PFClusterTools/interface/PFEnergyCalibration.h"
13 #include "RecoParticleFlow/PFClusterTools/interface/PFPhotonClusters.h"
14 #include "RecoParticleFlow/PFClusterTools/interface/PFEnergyCalibration.h"
15 #include "RecoParticleFlow/PFClusterTools/interface/PFSCEnergyCalibration.h"
16 #include "RecoParticleFlow/PFClusterTools/interface/PFEnergyResolution.h"
17 #include "RecoParticleFlow/PFClusterTools/interface/PFClusterWidthAlgo.h"
18 #include "RecoParticleFlow/PFProducer/interface/PFElectronExtraEqual.h"
19 #include "DataFormats/Common/interface/RefToPtr.h"
20 #include "RecoEcal/EgammaCoreTools/interface/Mustache.h"
21 #include "DataFormats/Math/interface/deltaPhi.h"
22 #include "DataFormats/Math/interface/deltaR.h"
23 #include <TFile.h>
24 #include <TVector2.h>
25 #include <iomanip>
26 #include <algorithm>
27 #include <TMath.h>
28 
29 // include combinations header (not yet included 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 
52 namespace {
53  typedef PFEGammaAlgo::PFSCElement SCElement;
54  typedef PFEGammaAlgo::EEtoPSAssociation EEtoPSAssociation;
55  typedef std::pair<CaloClusterPtr::key_type,CaloClusterPtr> EEtoPSElement;
56  typedef PFEGammaAlgo::PFClusterElement ClusterElement;
57  typedef PFEGammaAlgo::PFFlaggedElement PFFlaggedElement;
58  typedef PFEGammaAlgo::PFSCFlaggedElement SCFlaggedElement;
59  typedef PFEGammaAlgo::PFKFFlaggedElement KFFlaggedElement;
60  typedef PFEGammaAlgo::PFGSFFlaggedElement GSFFlaggedElement;
61  typedef PFEGammaAlgo::PFClusterFlaggedElement ClusterFlaggedElement;
62 
63  typedef std::unary_function<const ClusterFlaggedElement&,
64  bool> ClusterMatcher;
65 
66  typedef std::unary_function<const PFFlaggedElement&,
67  bool> PFFlaggedElementMatcher;
68  typedef std::binary_function<const PFFlaggedElement&,
69  const PFFlaggedElement&,
70  bool> PFFlaggedElementSorter;
71 
72  typedef std::unary_function<const reco::PFBlockElement&,
73  bool> PFElementMatcher;
74 
75  typedef std::unary_function<const PFEGammaAlgo::ProtoEGObject&,
76  bool> POMatcher;
77 
78  typedef std::unary_function<PFFlaggedElement&,
79  ClusterFlaggedElement> ClusterElementConverter;
80 
81  struct SumPSEnergy : public std::binary_function<double,
82  const ClusterFlaggedElement&,
83  double> {
84  reco::PFBlockElement::Type _thetype;
85  SumPSEnergy(reco::PFBlockElement::Type type) : _thetype(type) {}
86  double operator()(double a,
87  const ClusterFlaggedElement& b) {
88 
89  return a + (_thetype == b.first->type())*b.first->clusterRef()->energy();
90  }
91  };
92 
93  bool comparePSMapByKey(const EEtoPSElement& a,
94  const EEtoPSElement& b) {
95  return a.first < b.first;
96  }
97 
98  struct UsableElementToPSCluster : public ClusterElementConverter {
99  ClusterFlaggedElement operator () (PFFlaggedElement& elem) {
100  const ClusterElement* pselemascluster =
101  docast(const ClusterElement*,elem.first);
102  if( reco::PFBlockElement::PS1 != pselemascluster->type() &&
103  reco::PFBlockElement::PS2 != pselemascluster->type() ) {
104  std::stringstream ps_err;
105  pselemascluster->Dump(ps_err,"\t");
106  throw cms::Exception("UseableElementToPSCluster()")
107  << "This element is not a PS cluster!" << std::endl
108  << ps_err.str() << std::endl;
109  }
110  if( elem.second == false ) {
111  std::stringstream ps_err;
112  pselemascluster->Dump(ps_err,"\t");
113  throw cms::Exception("UsableElementToPSCluster()")
114  << "PS Cluster matched to EE is already used! "
115  << "This should be impossible!" << std::endl
116  << ps_err.str() << std::endl;
117  }
118  elem.second = false; // flag as used!
119  return std::make_pair(pselemascluster,true);
120  }
121  };
122 
123  struct SeedMatchesToSCElement : public PFFlaggedElementMatcher {
124  reco::SuperClusterRef _scfromseed;
125  SeedMatchesToSCElement(const reco::ElectronSeedRef& s) {
126  _scfromseed = s->caloCluster().castTo<reco::SuperClusterRef>();
127  }
128  bool operator() (const PFFlaggedElement& elem) {
129  const SCElement* scelem = docast(const SCElement*,elem.first);
130  if( _scfromseed.isNull() || !elem.second || !scelem) return false;
131  return ( _scfromseed->seed()->seed() ==
132  scelem->superClusterRef()->seed()->seed() );
133  }
134  };
135 
136  struct SCSubClusterMatchesToElement : public PFFlaggedElementMatcher {
137  const reco::CaloCluster_iterator begin, end;
138  SCSubClusterMatchesToElement(const reco::CaloCluster_iterator& b,
139  const reco::CaloCluster_iterator& e) :
140  begin(b),
141  end(e) { }
142  bool operator() (const PFFlaggedElement& elem) {
143  const ClusterElement* cluselem =
144  docast(const ClusterElement*,elem.first);
145  if( !elem.second || !cluselem) return false;
146  reco::CaloCluster_iterator cl = begin;
147  for( ; cl != end; ++cl ) {
148  if((*cl)->seed() == cluselem->clusterRef()->seed()) {
149  return true;
150  }
151  }
152  return false;
153  }
154  };
155 
156  struct SeedMatchesToProtoObject : public POMatcher {
157  reco::SuperClusterRef _scfromseed;
158  bool _ispfsc;
159  SeedMatchesToProtoObject(const reco::ElectronSeedRef& s) {
160  _scfromseed = s->caloCluster().castTo<reco::SuperClusterRef>();
161  _ispfsc = false;
162  if( _scfromseed.isNonnull() ) {
163  const edm::Ptr<reco::PFCluster> testCast(_scfromseed->seed());
164  _ispfsc = testCast.isNonnull();
165  }
166  }
167  bool operator() (const PFEGammaAlgo::ProtoEGObject& po) {
168  if( _scfromseed.isNull() || !po.parentSC ) return false;
169  if( _ispfsc ) {
170  return ( _scfromseed->seed() ==
171  po.parentSC->superClusterRef()->seed() );
172  }
173  return ( _scfromseed->seed()->seed() ==
174  po.parentSC->superClusterRef()->seed()->seed() );
175  }
176  };
177 
178  template<bool useConvs=false>
179  bool elementNotCloserToOther(const reco::PFBlockRef& block,
180  const PFBlockElement::Type& keytype,
181  const size_t key,
182  const PFBlockElement::Type& valtype,
183  const size_t test,
184  const float EoPin_cut = 1.0e6) {
185  constexpr reco::PFBlockElement::TrackType ConvType =
186  reco::PFBlockElement::T_FROM_GAMMACONV;
187  // this is inside out but I just want something that works right now
188  switch( keytype ) {
189  case reco::PFBlockElement::GSF:
190  {
191  const reco::PFBlockElementGsfTrack* elemasgsf =
192  docast(const reco::PFBlockElementGsfTrack*,
193  &(block->elements()[key]));
194  if( elemasgsf && valtype == PFBlockElement::ECAL ) {
195  const ClusterElement* elemasclus =
196  reinterpret_cast<const ClusterElement*>(&(block->elements()[test]));
197  float cluster_e = elemasclus->clusterRef()->correctedEnergy();
198  float trk_pin = elemasgsf->Pin().P();
199  if( cluster_e / trk_pin > EoPin_cut ) {
200  LOGDRESSED("elementNotCloserToOther")
201  << "GSF track failed EoP cut to match with cluster!";
202  return false;
203  }
204  }
205  }
206  break;
207  case reco::PFBlockElement::TRACK:
208  {
209  const reco::PFBlockElementTrack* elemaskf =
210  docast(const reco::PFBlockElementTrack*,
211  &(block->elements()[key]));
212  if( elemaskf && valtype == PFBlockElement::ECAL ) {
213  const ClusterElement* elemasclus =
214  reinterpret_cast<const ClusterElement*>(&(block->elements()[test]));
215  float cluster_e = elemasclus->clusterRef()->correctedEnergy();
216  float trk_pin =
217  std::sqrt(elemaskf->trackRef()->innerMomentum().mag2());
218  if( cluster_e / trk_pin > EoPin_cut ) {
219  LOGDRESSED("elementNotCloserToOther")
220  << "KF track failed EoP cut to match with cluster!";
221  return false;
222  }
223  }
224  }
225  break;
226  default:
227  break;
228  }
229 
230  const float dist =
231  block->dist(key,test,block->linkData(),reco::PFBlock::LINKTEST_ALL);
232  if( dist == -1.0f ) return false; // don't associate non-linked elems
233  std::multimap<double, unsigned> dists_to_val;
234  block->associatedElements(test,block->linkData(),dists_to_val,keytype,
235  reco::PFBlock::LINKTEST_ALL);
236  for( const auto& valdist : dists_to_val ) {
237  const size_t idx = valdist.second;
238  // check track types for conversion info
239  switch( keytype ) {
240  case reco::PFBlockElement::GSF:
241  {
242  const reco::PFBlockElementGsfTrack* elemasgsf =
243  docast(const reco::PFBlockElementGsfTrack*,
244  &(block->elements()[idx]));
245  if( !useConvs && elemasgsf->trackType(ConvType) ) return false;
246  if( elemasgsf && valtype == PFBlockElement::ECAL ) {
247  const ClusterElement* elemasclus =
248  docast(const ClusterElement*,&(block->elements()[test]));
249  float cluster_e = elemasclus->clusterRef()->correctedEnergy();
250  float trk_pin = elemasgsf->Pin().P();
251  if( cluster_e / trk_pin > EoPin_cut ) continue;
252  }
253  }
254  break;
255  case reco::PFBlockElement::TRACK:
256  {
257  const reco::PFBlockElementTrack* elemaskf =
258  docast(const reco::PFBlockElementTrack*,
259  &(block->elements()[idx]));
260  if( !useConvs && elemaskf->trackType(ConvType) ) return false;
261  if( elemaskf && valtype == PFBlockElement::ECAL ) {
262  const ClusterElement* elemasclus =
263  reinterpret_cast<const ClusterElement*>(&(block->elements()[test]));
264  float cluster_e = elemasclus->clusterRef()->correctedEnergy();
265  float trk_pin =
266  std::sqrt(elemaskf->trackRef()->innerMomentum().mag2());
267  if( cluster_e / trk_pin > EoPin_cut ) continue;
268  }
269  }
270  break;
271  default:
272  break;
273  }
274  if( valdist.first < dist && idx != key ) {
275  LOGDRESSED("elementNotCloserToOther")
276  << "key element of type " << keytype
277  << " is closer to another element of type" << valtype
278  << std::endl;
279  return false; // false if closer element of specified type found
280  }
281  }
282  return true;
283  }
284 
285  struct CompatibleEoPOut : public PFFlaggedElementMatcher {
286  const reco::PFBlockElementGsfTrack* comp;
287  CompatibleEoPOut(const reco::PFBlockElementGsfTrack* c) : comp(c) {}
288  bool operator()(const PFFlaggedElement& e) const {
289  if( PFBlockElement::ECAL != e.first->type() ) return false;
290  const ClusterElement* elemascluster =
291  docast(const ClusterElement*,e.first);
292  const float gsf_eta_diff = std::abs(comp->positionAtECALEntrance().eta()-
293  comp->Pout().eta());
294  const reco::PFClusterRef cRef = elemascluster->clusterRef();
295  return ( gsf_eta_diff <= 0.3 && cRef->energy()/comp->Pout().t() <= 5 );
296  }
297  };
298 
299  template<class TrackElementType>
300  struct IsConversionTrack : PFFlaggedElementMatcher {
301  bool operator()(const PFFlaggedElement& e) {
302  constexpr reco::PFBlockElement::TrackType ConvType =
303  reco::PFBlockElement::T_FROM_GAMMACONV;
304  const TrackElementType* elemastrk =
305  docast(const TrackElementType*,e.first);
306  return elemastrk->trackType(ConvType);
307  }
308  };
309 
310  template<PFBlockElement::Type keytype,
311  PFBlockElement::Type valtype,
312  bool useConv=false>
313  struct NotCloserToOther : public PFFlaggedElementMatcher {
314  const reco::PFBlockElement* comp;
315  const reco::PFBlockRef& block;
316  const reco::PFBlock::LinkData& links;
317  const float EoPin_cut;
318  NotCloserToOther(const reco::PFBlockRef& b,
319  const reco::PFBlock::LinkData& l,
320  const PFFlaggedElement* e,
321  const float EoPcut=1.0e6): comp(e->first),
322  block(b),
323  links(l),
324  EoPin_cut(EoPcut) {
325  }
326  NotCloserToOther(const reco::PFBlockRef& b,
327  const reco::PFBlock::LinkData& l,
328  const reco::PFBlockElement* e,
329  const float EoPcut=1.0e6): comp(e),
330  block(b),
331  links(l),
332  EoPin_cut(EoPcut) {
333  }
334  bool operator () (const PFFlaggedElement& e) {
335  if( !e.second || valtype != e.first->type() ) return false;
336  return elementNotCloserToOther<useConv>(block,
337  keytype,comp->index(),
338  valtype,e.first->index(),
339  EoPin_cut);
340  }
341  };
342 
343  struct LesserByDistance : public PFFlaggedElementSorter {
344  const reco::PFBlockElement* comp;
345  const reco::PFBlockRef& block;
346  const reco::PFBlock::LinkData& links;
347  LesserByDistance(const reco::PFBlockRef& b,
348  const reco::PFBlock::LinkData& l,
349  const PFFlaggedElement* e): comp(e->first),
350  block(b),
351  links(l) {}
352  LesserByDistance(const reco::PFBlockRef& b,
353  const reco::PFBlock::LinkData& l,
354  const reco::PFBlockElement* e): comp(e),
355  block(b),
356  links(l) {}
357  bool operator () (const PFFlaggedElement& e1,
358  const PFFlaggedElement& e2) {
359  double dist1 = block->dist(comp->index(),
360  e1.first->index(),
361  links,
362  reco::PFBlock::LINKTEST_ALL);
363  double dist2 = block->dist(comp->index(),
364  e2.first->index(),
365  links,
366  reco::PFBlock::LINKTEST_ALL);
367  dist1 = ( dist1 == -1.0 ? 1e6 : dist1 );
368  dist2 = ( dist2 == -1.0 ? 1e6 : dist2 );
369  return dist1 < dist2;
370  }
371  };
372 
373  bool isROLinkedByClusterOrTrack(const PFEGammaAlgo::ProtoEGObject& RO1,
374  const PFEGammaAlgo::ProtoEGObject& RO2 ) {
375  // also don't allow ROs where both have clusters
376  // and GSF tracks to merge (10 Dec 2013)
377  if(RO1.primaryGSFs.size() && RO2.primaryGSFs.size()) {
378  LOGDRESSED("isROLinkedByClusterOrTrack")
379  << "cannot merge, both have GSFs!" << std::endl;
380  return false;
381  }
382  // don't allow EB/EE to mix (11 Sept 2013)
383  if( RO1.ecalclusters.size() && RO2.ecalclusters.size() ) {
384  if(RO1.ecalclusters.front().first->clusterRef()->layer() !=
385  RO2.ecalclusters.front().first->clusterRef()->layer() ) {
386  LOGDRESSED("isROLinkedByClusterOrTrack")
387  << "cannot merge, different ECAL types!" << std::endl;
388  return false;
389  }
390  }
391  const reco::PFBlockRef& blk = RO1.parentBlock;
392  bool not_closer;
393  // check links track -> cluster
394  for( const auto& cluster: RO1.ecalclusters ) {
395  for( const auto& primgsf : RO2.primaryGSFs ) {
396  not_closer =
397  elementNotCloserToOther(blk,
398  cluster.first->type(),
399  cluster.first->index(),
400  primgsf.first->type(),
401  primgsf.first->index());
402  if( not_closer ) {
403  LOGDRESSED("isROLinkedByClusterOrTrack")
404  << "merged by cluster to primary GSF" << std::endl;
405  return true;
406  } else {
407  LOGDRESSED("isROLinkedByClusterOrTrack")
408  << "cluster to primary GSF failed since"
409  << " cluster closer to another GSF" << std::endl;
410  }
411  }
412  for( const auto& primkf : RO2.primaryKFs) {
413  not_closer =
414  elementNotCloserToOther(blk,
415  cluster.first->type(),
416  cluster.first->index(),
417  primkf.first->type(),
418  primkf.first->index());
419  if( not_closer ) {
420  LOGDRESSED("isROLinkedByClusterOrTrack")
421  << "merged by cluster to primary KF" << std::endl;
422  return true;
423  }
424  }
425  for( const auto& secdkf : RO2.secondaryKFs) {
426  not_closer =
427  elementNotCloserToOther(blk,
428  cluster.first->type(),
429  cluster.first->index(),
430  secdkf.first->type(),
431  secdkf.first->index());
432  if( not_closer ) {
433  LOGDRESSED("isROLinkedByClusterOrTrack")
434  << "merged by cluster to secondary KF" << std::endl;
435  return true;
436  }
437  }
438  // check links brem -> cluster
439  for( const auto& brem : RO2.brems ) {
440  not_closer = elementNotCloserToOther(blk,
441  cluster.first->type(),
442  cluster.first->index(),
443  brem.first->type(),
444  brem.first->index());
445  if( not_closer ) {
446  LOGDRESSED("isROLinkedByClusterOrTrack")
447  << "merged by cluster to brem KF" << std::endl;
448  return true;
449  }
450  }
451  }
452  // check links primary gsf -> secondary kf
453  for( const auto& primgsf : RO1.primaryGSFs ) {
454  for( const auto& secdkf : RO2.secondaryKFs) {
455  not_closer =
456  elementNotCloserToOther(blk,
457  primgsf.first->type(),
458  primgsf.first->index(),
459  secdkf.first->type(),
460  secdkf.first->index());
461  if( not_closer ) {
462  LOGDRESSED("isROLinkedByClusterOrTrack")
463  << "merged by GSF to secondary KF" << std::endl;
464  return true;
465  }
466  }
467  }
468  // check links primary kf -> secondary kf
469  for( const auto& primkf : RO1.primaryKFs ) {
470  for( const auto& secdkf : RO2.secondaryKFs) {
471  not_closer =
472  elementNotCloserToOther(blk,
473  primkf.first->type(),
474  primkf.first->index(),
475  secdkf.first->type(),
476  secdkf.first->index());
477  if( not_closer ) {
478  LOGDRESSED("isROLinkedByClusterOrTrack")
479  << "merged by primary KF to secondary KF" << std::endl;
480  return true;
481  }
482  }
483  }
484  // check links secondary kf -> secondary kf
485  for( const auto& secdkf1 : RO1.secondaryKFs ) {
486  for( const auto& secdkf2 : RO2.secondaryKFs) {
487  not_closer =
488  elementNotCloserToOther<true>(blk,
489  secdkf1.first->type(),
490  secdkf1.first->index(),
491  secdkf2.first->type(),
492  secdkf2.first->index());
493  if( not_closer ) {
494  LOGDRESSED("isROLinkedByClusterOrTrack")
495  << "merged by secondary KF to secondary KF" << std::endl;
496  return true;
497  }
498  }
499  }
500  return false;
501  }
502 
503  struct TestIfROMergableByLink : public POMatcher {
504  const PFEGammaAlgo::ProtoEGObject& comp;
505  TestIfROMergableByLink(const PFEGammaAlgo::ProtoEGObject& RO) :
506  comp(RO) {}
507  bool operator() (const PFEGammaAlgo::ProtoEGObject& ro) {
508  const bool result = ( isROLinkedByClusterOrTrack(comp,ro) ||
509  isROLinkedByClusterOrTrack(ro,comp) );
510  return result;
511  }
512  };
513 
514  std::vector<const ClusterElement*>
515  getSCAssociatedECALsSafe(const reco::SuperClusterRef& scref,
516  std::vector<PFFlaggedElement>& ecals) {
517  std::vector<const ClusterElement*> cluster_list;
518  auto sccl = scref->clustersBegin();
519  auto scend = scref->clustersEnd();
520  auto pfc = ecals.begin();
521  auto pfcend = ecals.end();
522  for( ; sccl != scend; ++sccl ) {
523  std::vector<const ClusterElement*> matched_pfcs;
524  const double eg_energy = (*sccl)->energy();
525 
526  for( pfc = ecals.begin(); pfc != pfcend; ++pfc ) {
527  const ClusterElement *pfcel =
528  docast(const ClusterElement*, pfc->first);
529  const bool matched =
530  ClusterClusterMapping::overlap(**sccl,*(pfcel->clusterRef()));
531  // need to protect against high energy clusters being attached
532  // to low-energy SCs
533  if( matched && pfcel->clusterRef()->energy() < 1.2*scref->energy()) {
534  matched_pfcs.push_back(pfcel);
535  }
536  }
537  std::sort(matched_pfcs.begin(),matched_pfcs.end());
538 
539  double min_residual = 1e6;
540  std::vector<const ClusterElement*> best_comb;
541  for( size_t i = 1; i <= matched_pfcs.size(); ++i ) {
542  //now we find the combination of PF-clusters which
543  //has the smallest energy residual with respect to the
544  //EG-cluster we are looking at now
545  do {
546  double energy = std::accumulate(matched_pfcs.begin(),
547  matched_pfcs.begin()+i,
548  0.0,
549  [](const double a,
550  const ClusterElement* c)
551  { return a + c->clusterRef()->energy(); });
552  const double resid = std::abs(energy - eg_energy);
553  if( resid < min_residual ) {
554  best_comb.clear();
555  best_comb.reserve(i);
556  min_residual = resid;
557  best_comb.insert(best_comb.begin(),
558  matched_pfcs.begin(),
559  matched_pfcs.begin()+i);
560  }
561  }while(boost::next_combination(matched_pfcs.begin(),
562  matched_pfcs.begin()+i,
563  matched_pfcs.end()));
564  }
565  for( const auto& clelem : best_comb ) {
566  reco::PFClusterRef clref = clelem->clusterRef();
567  if( std::find(cluster_list.begin(),cluster_list.end(),clelem) ==
568  cluster_list.end() ) {
569  cluster_list.push_back(clelem);
570  }
571  }
572  }
573  return cluster_list;
574  }
575  bool addPFClusterToROSafe(const ClusterElement* cl,
576  PFEGammaAlgo::ProtoEGObject& RO) {
577  if( !RO.ecalclusters.size() ) {
578  RO.ecalclusters.push_back(std::make_pair(cl,true));
579  return true;
580  } else {
581  const PFLayer::Layer clayer = cl->clusterRef()->layer();
582  const PFLayer::Layer blayer =
583  RO.ecalclusters.back().first->clusterRef()->layer();
584  if( clayer == blayer ) {
585  RO.ecalclusters.push_back(std::make_pair(cl,true));
586  return true;
587  }
588  }
589  return false;
590  }
591 
592  // sets the cluster best associated to the GSF track
593  // leave it null if no GSF track
594  void setROElectronCluster(PFEGammaAlgo::ProtoEGObject& RO) {
595  if( !RO.ecalclusters.size() ) return;
596  RO.lateBrem = -1;
597  RO.firstBrem = -1;
598  RO.nBremsWithClusters = -1;
599  const reco::PFBlockElementBrem *firstBrem = NULL, *lastBrem = NULL;
600  const reco::PFBlockElementCluster *bremCluster = NULL, *gsfCluster = NULL,
601  *kfCluster = NULL, *gsfCluster_noassc = NULL;
602  const reco::PFBlockRef& parent = RO.parentBlock;
603  int nBremClusters = 0;
604  constexpr float maxDist = 1e6;
605  float mDist_gsf(maxDist), mDist_gsf_noassc(maxDist), mDist_kf(maxDist);
606  for( const auto& cluster : RO.ecalclusters ) {
607  for( const auto& gsf : RO.primaryGSFs ) {
608  const bool hasclu = elementNotCloserToOther(parent,
609  gsf.first->type(),
610  gsf.first->index(),
611  cluster.first->type(),
612  cluster.first->index());
613  const float deta =
614  std::abs(cluster.first->clusterRef()->positionREP().eta() -
615  gsf.first->positionAtECALEntrance().eta());
616  const float dphi =
617  std::abs(TVector2::Phi_mpi_pi(
618  cluster.first->clusterRef()->positionREP().phi() -
619  gsf.first->positionAtECALEntrance().phi()));
620  const float dist = std::hypot(deta,dphi);
621  if( hasclu && dist < mDist_gsf ) {
622  gsfCluster = cluster.first;
623  mDist_gsf = dist;
624  } else if ( dist < mDist_gsf_noassc ) {
625  gsfCluster_noassc = cluster.first;
626  mDist_gsf_noassc = dist;
627  }
628  }
629  for( const auto& kf : RO.primaryKFs ) {
630  const bool hasclu = elementNotCloserToOther(parent,
631  kf.first->type(),
632  kf.first->index(),
633  cluster.first->type(),
634  cluster.first->index());
635  const float dist = parent->dist(cluster.first->index(),
636  kf.first->index(),
637  parent->linkData(),
638  reco::PFBlock::LINKTEST_ALL);
639  if( hasclu && dist < mDist_kf ) {
640  kfCluster = cluster.first;
641  mDist_kf = dist;
642  }
643  }
644  for( const auto& brem : RO.brems ) {
645  const bool hasclu = elementNotCloserToOther(parent,
646  brem.first->type(),
647  brem.first->index(),
648  cluster.first->type(),
649  cluster.first->index());
650  if( hasclu ) {
651  ++nBremClusters;
652  if( !firstBrem ||
653  ( firstBrem->indTrajPoint() - 2 >
654  brem.first->indTrajPoint() - 2) ) {
655  firstBrem = brem.first;
656  }
657  if( !lastBrem ||
658  ( lastBrem->indTrajPoint() - 2 <
659  brem.first->indTrajPoint() - 2) ) {
660  lastBrem = brem.first;
661  bremCluster = cluster.first;
662  }
663  }
664  }
665  }
666  if( !gsfCluster && !kfCluster && !bremCluster ) {
667  gsfCluster = gsfCluster_noassc;
668  }
669  RO.nBremsWithClusters = nBremClusters;
670  RO.lateBrem = 0;
671  if( gsfCluster ) {
672  RO.electronClusters.push_back(gsfCluster);
673  } else if ( kfCluster ) {
674  RO.electronClusters.push_back(kfCluster);
675  }
676  if( bremCluster && !gsfCluster && !kfCluster ) {
677  RO.electronClusters.push_back(bremCluster);
678  }
679  if( firstBrem && RO.ecalclusters.size() > 1 ) {
680  RO.firstBrem = firstBrem->indTrajPoint() - 2;
681  if( bremCluster == gsfCluster ) RO.lateBrem = 1;
682  }
683  }
684 }
685 
686 PFEGammaAlgo::
687 PFEGammaAlgo(const PFEGammaAlgo::PFEGConfigInfo& cfg) :
688  cfg_(cfg),
689  isvalid_(false),
690  verbosityLevel_(Silent),
691  nlost(0.0), nlayers(0.0),
692  chi2(0.0), STIP(0.0), del_phi(0.0),HoverPt(0.0), EoverPt(0.0), track_pt(0.0),
693  mvaValue(0.0),
694  CrysPhi_(0.0), CrysEta_(0.0), VtxZ_(0.0), ClusPhi_(0.0), ClusEta_(0.0),
695  ClusR9_(0.0), Clus5x5ratio_(0.0), PFCrysEtaCrack_(0.0), logPFClusE_(0.0), e3x3_(0.0),
696  CrysIPhi_(0), CrysIEta_(0),
697  CrysX_(0.0), CrysY_(0.0),
698  EB(0.0),
699  eSeed_(0.0), e1x3_(0.0),e3x1_(0.0), e1x5_(0.0), e2x5Top_(0.0), e2x5Bottom_(0.0), e2x5Left_(0.0), e2x5Right_(0.0),
700  etop_(0.0), ebottom_(0.0), eleft_(0.0), eright_(0.0),
701  e2x5Max_(0.0),
702  PFPhoEta_(0.0), PFPhoPhi_(0.0), PFPhoR9_(0.0), PFPhoR9Corr_(0.0), SCPhiWidth_(0.0), SCEtaWidth_(0.0),
703  PFPhoEt_(0.0), RConv_(0.0), PFPhoEtCorr_(0.0), PFPhoE_(0.0), PFPhoECorr_(0.0), MustE_(0.0), E3x3_(0.0),
704  dEta_(0.0), dPhi_(0.0), LowClusE_(0.0), RMSAll_(0.0), RMSMust_(0.0), nPFClus_(0.0),
705  TotPS1_(0.0), TotPS2_(0.0),
706  nVtx_(0.0),
707  x0inner_(0.0), x0middle_(0.0), x0outer_(0.0),
708  excluded_(0.0), Mustache_EtRatio_(0.0), Mustache_Et_out_(0.0)
709 {
710 
711  // Set the tmva reader for electrons
712  tmvaReaderEle_ = new TMVA::Reader("!Color:Silent");
713  tmvaReaderEle_->AddVariable("lnPt_gsf",&lnPt_gsf);
714  tmvaReaderEle_->AddVariable("Eta_gsf",&Eta_gsf);
715  tmvaReaderEle_->AddVariable("dPtOverPt_gsf",&dPtOverPt_gsf);
716  tmvaReaderEle_->AddVariable("DPtOverPt_gsf",&DPtOverPt_gsf);
717  //tmvaReaderEle_->AddVariable("nhit_gsf",&nhit_gsf);
718  tmvaReaderEle_->AddVariable("chi2_gsf",&chi2_gsf);
719  //tmvaReaderEle_->AddVariable("DPtOverPt_kf",&DPtOverPt_kf);
720  tmvaReaderEle_->AddVariable("nhit_kf",&nhit_kf);
721  tmvaReaderEle_->AddVariable("chi2_kf",&chi2_kf);
722  tmvaReaderEle_->AddVariable("EtotPinMode",&EtotPinMode);
723  tmvaReaderEle_->AddVariable("EGsfPoutMode",&EGsfPoutMode);
724  tmvaReaderEle_->AddVariable("EtotBremPinPoutMode",&EtotBremPinPoutMode);
725  tmvaReaderEle_->AddVariable("DEtaGsfEcalClust",&DEtaGsfEcalClust);
726  tmvaReaderEle_->AddVariable("SigmaEtaEta",&SigmaEtaEta);
727  tmvaReaderEle_->AddVariable("HOverHE",&HOverHE);
728 // tmvaReaderEle_->AddVariable("HOverPin",&HOverPin);
729  tmvaReaderEle_->AddVariable("lateBrem",&lateBrem);
730  tmvaReaderEle_->AddVariable("firstBrem",&firstBrem);
731  tmvaReaderEle_->BookMVA("BDT",cfg_.mvaWeightFileEleID.c_str());
732 
733 
734  //Book MVA
735  tmvaReader_ = new TMVA::Reader("!Color:Silent");
736  tmvaReader_->AddVariable("del_phi",&del_phi);
737  tmvaReader_->AddVariable("nlayers", &nlayers);
738  tmvaReader_->AddVariable("chi2",&chi2);
739  tmvaReader_->AddVariable("EoverPt",&EoverPt);
740  tmvaReader_->AddVariable("HoverPt",&HoverPt);
741  tmvaReader_->AddVariable("track_pt", &track_pt);
742  tmvaReader_->AddVariable("STIP",&STIP);
743  tmvaReader_->AddVariable("nlost", &nlost);
744  tmvaReader_->BookMVA("BDT",cfg_.mvaweightfile.c_str());
745 
746  //Material Map
747  TFile *XO_File = new TFile(cfg_.X0_Map.c_str(),"READ");
748  X0_sum = (TH2D*)XO_File->Get("TrackerSum");
749  X0_inner = (TH2D*)XO_File->Get("Inner");
750  X0_middle = (TH2D*)XO_File->Get("Middle");
751  X0_outer = (TH2D*)XO_File->Get("Outer");
752 
753 }
754 
755 void PFEGammaAlgo::RunPFEG(const reco::PFBlockRef& blockRef,
756  std::vector<bool>& active) {
757 
758  fifthStepKfTrack_.clear();
759  convGsfTrack_.clear();
760 
761  egCandidate_.clear();
762  egExtra_.clear();
763 
764  // define how much is printed out for debugging.
765  // ... will be setable via CFG file parameter
766  verbosityLevel_ = Chatty; // Chatty mode.
767 
768  buildAndRefineEGObjects(blockRef);
769 }
770 
771 float PFEGammaAlgo::EvaluateSingleLegMVA(const reco::PFBlockRef& blockref,
772  const reco::Vertex& primaryvtx,
773  unsigned int track_index) {
774  const reco::PFBlock& block = *blockref;
775  const edm::OwnVector< reco::PFBlockElement >& elements = block.elements();
776  //use this to store linkdata in the associatedElements function below
777  PFBlock::LinkData linkData = block.linkData();
778  //calculate MVA Variables
779  chi2=elements[track_index].trackRef()->chi2()/elements[track_index].trackRef()->ndof();
780  nlost=elements[track_index].trackRef()->trackerExpectedHitsInner().numberOfLostHits();
781  nlayers=elements[track_index].trackRef()->hitPattern().trackerLayersWithMeasurement();
782  track_pt=elements[track_index].trackRef()->pt();
783  STIP=elements[track_index].trackRefPF()->STIP();
784 
785  float linked_e=0;
786  float linked_h=0;
787  std::multimap<double, unsigned int> ecalAssoTrack;
788  block.associatedElements( track_index,linkData,
789  ecalAssoTrack,
790  reco::PFBlockElement::ECAL,
791  reco::PFBlock::LINKTEST_ALL );
792  std::multimap<double, unsigned int> hcalAssoTrack;
793  block.associatedElements( track_index,linkData,
794  hcalAssoTrack,
795  reco::PFBlockElement::HCAL,
796  reco::PFBlock::LINKTEST_ALL );
797  if(ecalAssoTrack.size() > 0) {
798  for(std::multimap<double, unsigned int>::iterator itecal = ecalAssoTrack.begin();
799  itecal != ecalAssoTrack.end(); ++itecal) {
800  linked_e=linked_e+elements[itecal->second].clusterRef()->energy();
801  }
802  }
803  if(hcalAssoTrack.size() > 0) {
804  for(std::multimap<double, unsigned int>::iterator ithcal = hcalAssoTrack.begin();
805  ithcal != hcalAssoTrack.end(); ++ithcal) {
806  linked_h=linked_h+elements[ithcal->second].clusterRef()->energy();
807  }
808  }
809  EoverPt=linked_e/elements[track_index].trackRef()->pt();
810  HoverPt=linked_h/elements[track_index].trackRef()->pt();
811  GlobalVector rvtx(elements[track_index].trackRef()->innerPosition().X()-primaryvtx.x(),
812  elements[track_index].trackRef()->innerPosition().Y()-primaryvtx.y(),
813  elements[track_index].trackRef()->innerPosition().Z()-primaryvtx.z());
814  double vtx_phi=rvtx.phi();
815  //delta Phi between conversion vertex and track
816  del_phi=fabs(deltaPhi(vtx_phi, elements[track_index].trackRef()->innerMomentum().Phi()));
817  mvaValue = tmvaReader_->EvaluateMVA("BDT");
818 
819  return mvaValue;
820 }
821 
822 bool PFEGammaAlgo::isAMuon(const reco::PFBlockElement& pfbe) {
823  switch( pfbe.type() ) {
824  case reco::PFBlockElement::GSF:
825  {
826  auto& elements = _currentblock->elements();
827  std::multimap<double,unsigned> tks;
828  _currentblock->associatedElements(pfbe.index(),
829  _currentlinks,
830  tks,
831  reco::PFBlockElement::TRACK,
832  reco::PFBlock::LINKTEST_ALL);
833  for( const auto& tk : tks ) {
834  if( PFMuonAlgo::isMuon(elements[tk.second]) ) {
835  return true;
836  }
837  }
838  }
839  break;
840  case reco::PFBlockElement::TRACK:
841  return PFMuonAlgo::isMuon(pfbe);
842  break;
843  default:
844  break;
845  }
846  return false;
847 }
848 
849 void PFEGammaAlgo::buildAndRefineEGObjects(const reco::PFBlockRef& block) {
850  LOGVERB("PFEGammaAlgo")
851  << "Resetting PFEGammaAlgo for new block and running!" << std::endl;
852  _splayedblock.clear();
853  _recoveredlinks.clear();
854  _refinableObjects.clear();
855  _finalCandidates.clear();
856  _splayedblock.resize(12); // make sure that we always have the SC entry
857 
858  _currentblock = block;
859  _currentlinks = block->linkData();
860  //LOGDRESSED("PFEGammaAlgo") << *_currentblock << std::endl;
861  LOGVERB("PFEGammaAlgo") << "Splaying block" << std::endl;
862  //unwrap the PF block into a fast access map
863  for( const auto& pfelement : _currentblock->elements() ) {
864  if( isAMuon(pfelement) ) continue; // don't allow muons in our element list
865  const size_t itype = (size_t)pfelement.type();
866  if( itype >= _splayedblock.size() ) _splayedblock.resize(itype+1);
867  _splayedblock[itype].push_back(std::make_pair(&pfelement,true));
868  }
869 
870  // show the result of splaying the tree if it's really *really* needed
871 #ifdef PFLOW_DEBUG
872  std::stringstream splayout;
873  for( size_t itype = 0; itype < _splayedblock.size(); ++itype ) {
874  splayout << "\tType: " << itype << " indices: ";
875  for( const auto& flaggedelement : _splayedblock[itype] ) {
876  splayout << flaggedelement.first->index() << ' ';
877  }
878  if( itype != _splayedblock.size() - 1 ) splayout << std::endl;
879  }
880  LOGVERB("PFEGammaAlgo") << splayout.str();
881 #endif
882 
883  // precleaning of the ECAL clusters with respect to primary KF tracks
884  // we don't allow clusters in super clusters to be locked out this way
885  removeOrLinkECALClustersToKFTracks();
886 
887  initializeProtoCands(_refinableObjects);
888  LOGDRESSED("PFEGammaAlgo")
889  << "Initialized " << _refinableObjects.size() << " proto-EGamma objects"
890  << std::endl;
891  dumpCurrentRefinableObjects();
892 
893  //
894  // now we start the refining steps
895  //
896  //
897 
898  // --- Primary Linking Step ---
899  // since this is particle flow and we try to work from the pixels out
900  // we start by linking the tracks together and finding the ECAL clusters
901  for( auto& RO : _refinableObjects ) {
902  // find the KF tracks associated to GSF primary tracks
903  linkRefinableObjectGSFTracksToKFs(RO);
904  // do the same for HCAL clusters associated to the GSF
905  linkRefinableObjectPrimaryGSFTrackToHCAL(RO);
906  // link secondary KF tracks associated to primary KF tracks
907  linkRefinableObjectPrimaryKFsToSecondaryKFs(RO);
908  // pick up clusters that are linked to the GSF primary
909  linkRefinableObjectPrimaryGSFTrackToECAL(RO);
910  // link associated KF to ECAL (ECAL part grabs PS clusters too if able)
911  linkRefinableObjectKFTracksToECAL(RO);
912  // now finally look for clusters associated to brem tangents
913  linkRefinableObjectBremTangentsToECAL(RO);
914  }
915 
916  LOGDRESSED("PFEGammaAlgo")
917  << "Dumping after GSF and KF Track (Primary) Linking : " << std::endl;
918  dumpCurrentRefinableObjects();
919 
920  // merge objects after primary linking
921  mergeROsByAnyLink(_refinableObjects);
922 
923  LOGDRESSED("PFEGammaAlgo")
924  << "Dumping after first merging operation : " << std::endl;
925  dumpCurrentRefinableObjects();
926 
927  // --- Secondary Linking Step ---
928  // after this we go through the ECAL clusters on the remaining tracks
929  // and try to link those in...
930  for( auto& RO : _refinableObjects ) {
931  // look for conversion legs
932  linkRefinableObjectECALToSingleLegConv(RO);
933  dumpCurrentRefinableObjects();
934  // look for tracks that complement conversion legs
935  linkRefinableObjectConvSecondaryKFsToSecondaryKFs(RO);
936  // look again for ECAL clusters (this time with an e/p cut)
937  linkRefinableObjectSecondaryKFsToECAL(RO);
938  }
939 
940  LOGDRESSED("PFEGammaAlgo")
941  << "Dumping after ECAL to Track (Secondary) Linking : " << std::endl;
942  dumpCurrentRefinableObjects();
943 
944  // merge objects after primary linking
945  mergeROsByAnyLink(_refinableObjects);
946 
947  LOGDRESSED("PFEGammaAlgo")
948  << "There are " << _refinableObjects.size()
949  << " after the 2nd merging step." << std::endl;
950  dumpCurrentRefinableObjects();
951 
952  // -- unlinking and proto-object vetos, final sorting
953  for( auto& RO : _refinableObjects ) {
954  // remove secondary KFs (and possibly ECALs) matched to HCAL clusters
955  unlinkRefinableObjectKFandECALMatchedToHCAL(RO, false, false);
956  // remove secondary KFs and ECALs linked to them that have bad E/p_in
957  // and spoil the resolution
958  unlinkRefinableObjectKFandECALWithBadEoverP(RO);
959  // put things back in order after partitioning
960  std::sort(RO.ecalclusters.begin(), RO.ecalclusters.end(),
961  [](const PFClusterFlaggedElement& a,
962  const PFClusterFlaggedElement& b)
963  { return ( a.first->clusterRef()->correctedEnergy() >
964  b.first->clusterRef()->correctedEnergy() ) ; });
965  setROElectronCluster(RO);
966  }
967 
968  LOGDRESSED("PFEGammaAlgo")
969  << "There are " << _refinableObjects.size()
970  << " after the unlinking and vetos step." << std::endl;
971  dumpCurrentRefinableObjects();
972 
973  // fill the PF candidates and then build the refined SC
974  fillPFCandidates(_refinableObjects,outcands_,outcandsextra_);
975 
976 }
977 
978 void PFEGammaAlgo::
979 initializeProtoCands(std::list<PFEGammaAlgo::ProtoEGObject>& egobjs) {
980  // step 1: build SC based proto-candidates
981  // in the future there will be an SC Et requirement made here to control
982  // block size
983  for( auto& element : _splayedblock[PFBlockElement::SC] ) {
984  LOGDRESSED("PFEGammaAlgo")
985  << "creating SC-based proto-object" << std::endl
986  << "\tSC at index: " << element.first->index()
987  << " has type: " << element.first->type() << std::endl;
988  element.second = false;
989  ProtoEGObject fromSC;
990  fromSC.nBremsWithClusters = -1;
991  fromSC.firstBrem = -1;
992  fromSC.lateBrem = -1;
993  fromSC.parentBlock = _currentblock;
994  fromSC.parentSC = docast(const PFSCElement*,element.first);
995  // splay the supercluster so we can knock out used elements
996  bool sc_success =
997  unwrapSuperCluster(fromSC.parentSC,fromSC.ecalclusters,fromSC.ecal2ps);
998  if( sc_success ) {
999  /*
1000  auto ins_pos = std::lower_bound(_refinableObjects.begin(),
1001  _refinableObjects.end(),
1002  fromSC,
1003  [&](const ProtoEGObject& a,
1004  const ProtoEGObject& b){
1005  const double a_en =
1006  a.parentSC->superClusterRef()->energy();
1007  const double b_en =
1008  b.parentSC->superClusterRef()->energy();
1009  return a_en < b_en;
1010  });
1011  */
1012  _refinableObjects.insert(_refinableObjects.end(),fromSC);
1013  }
1014  }
1015  // step 2: build GSF-seed-based proto-candidates
1016  reco::GsfTrackRef gsfref_forextra;
1017  reco::TrackExtraRef gsftrk_extra;
1018  reco::ElectronSeedRef theseedref;
1019  std::list<ProtoEGObject>::iterator objsbegin, objsend;
1020  for( auto& element : _splayedblock[PFBlockElement::GSF] ) {
1021  LOGDRESSED("PFEGammaAlgo")
1022  << "creating GSF-based proto-object" << std::endl
1023  << "\tGSF at index: " << element.first->index()
1024  << " has type: " << element.first->type() << std::endl;
1025  const PFGSFElement* elementAsGSF =
1026  docast(const PFGSFElement*,element.first);
1027  if( elementAsGSF->trackType(reco::PFBlockElement::T_FROM_GAMMACONV) ) {
1028  continue; // for now, do not allow dedicated brems to make proto-objects
1029  }
1030  element.second = false;
1031 
1032  ProtoEGObject fromGSF;
1033  fromGSF.nBremsWithClusters = -1;
1034  fromGSF.firstBrem = -1;
1035  fromGSF.lateBrem = 0;
1036  gsfref_forextra = elementAsGSF->GsftrackRef();
1037  gsftrk_extra = ( gsfref_forextra.isAvailable() ?
1038  gsfref_forextra->extra() : reco::TrackExtraRef() );
1039  theseedref = ( gsftrk_extra.isAvailable() ?
1040  gsftrk_extra->seedRef().castTo<reco::ElectronSeedRef>() :
1041  reco::ElectronSeedRef() );
1042  fromGSF.electronSeed = theseedref;
1043  // exception if there's no seed
1044  if(fromGSF.electronSeed.isNull() || !fromGSF.electronSeed.isAvailable()) {
1045  std::stringstream gsf_err;
1046  elementAsGSF->Dump(gsf_err,"\t");
1047  throw cms::Exception("PFEGammaAlgo::initializeProtoCands()")
1048  << "Found a GSF track with no seed! This should not happen!"
1049  << std::endl << gsf_err.str() << std::endl;
1050  }
1051  // flag this GSF element as globally used and push back the track ref
1052  // into the protocand
1053  element.second = false;
1054  fromGSF.parentBlock = _currentblock;
1055  fromGSF.primaryGSFs.push_back(std::make_pair(elementAsGSF,true));
1056  // add the directly matched brem tangents
1057  for( auto& brem : _splayedblock[PFBlockElement::BREM] ) {
1058  float dist = _currentblock->dist(elementAsGSF->index(),
1059  brem.first->index(),
1060  _currentlinks,
1061  reco::PFBlock::LINKTEST_ALL);
1062  if( dist == 0.001f ) {
1063  const PFBremElement* eAsBrem =
1064  docast(const PFBremElement*,brem.first);
1065  fromGSF.brems.push_back(std::make_pair(eAsBrem,true));
1066  fromGSF.localMap.push_back( ElementMap::value_type(eAsBrem,elementAsGSF) );
1067  fromGSF.localMap.push_back( ElementMap::value_type(elementAsGSF,eAsBrem) );
1068  brem.second = false;
1069  }
1070  }
1071  // if this track is ECAL seeded reset links or import cluster
1072  // tracker (this is pixel only, right?) driven seeds just get the GSF
1073  // track associated since this only branches for ECAL Driven seeds
1074  if( fromGSF.electronSeed->isEcalDriven() ) {
1075  // step 2a: either merge with existing ProtoEG object with SC or add
1076  // SC directly to this proto EG object if not present
1077  LOGDRESSED("PFEGammaAlgo")
1078  << "GSF-based proto-object is ECAL driven, merging SC-cand"
1079  << std::endl;
1080  LOGVERB("PFEGammaAlgo")
1081  << "ECAL Seed Ptr: " << fromGSF.electronSeed.get()
1082  << " isAvailable: " << fromGSF.electronSeed.isAvailable()
1083  << " isNonnull: " << fromGSF.electronSeed.isNonnull()
1084  << std::endl;
1085  SeedMatchesToProtoObject sctoseedmatch(fromGSF.electronSeed);
1086  objsbegin = _refinableObjects.begin();
1087  objsend = _refinableObjects.end();
1088  // this auto is a std::list<ProtoEGObject>::iterator
1089  auto clusmatch = std::find_if(objsbegin,objsend,sctoseedmatch);
1090  if( clusmatch != objsend ) {
1091  fromGSF.parentSC = clusmatch->parentSC;
1092  fromGSF.ecalclusters = std::move(clusmatch->ecalclusters);
1093  fromGSF.ecal2ps = std::move(clusmatch->ecal2ps);
1094  _refinableObjects.erase(clusmatch);
1095  } else if (fromGSF.electronSeed.isAvailable() &&
1096  fromGSF.electronSeed.isNonnull()) {
1097  // link tests in the gap region can current split a gap electron
1098  // HEY THIS IS A WORK AROUND FOR A KNOWN BUG IN PFBLOCKALGO
1099  // MAYBE WE SHOULD FIX IT??????????????????????????????????
1100  LOGDRESSED("PFEGammaAlgo")
1101  << "Encountered the known GSF-SC splitting bug "
1102  << " in PFBlockAlgo! We should really fix this!" << std::endl;
1103  } else { // SC was not in a earlier proto-object
1104  std::stringstream gsf_err;
1105  elementAsGSF->Dump(gsf_err,"\t");
1106  throw cms::Exception("PFEGammaAlgo::initializeProtoCands()")
1107  << "Expected SuperCluster from ECAL driven GSF seed "
1108  << "was not found in the block!" << std::endl
1109  << gsf_err.str() << std::endl;
1110  } // supercluster in block
1111  } // is ECAL driven seed?
1112  /*
1113  auto ins_pos = std::lower_bound(_refinableObjects.begin(),
1114  _refinableObjects.end(),
1115  fromGSF,
1116  [&](const ProtoEGObject& a,
1117  const ProtoEGObject& b){
1118  const double a_en = ( a.parentSC ?
1119  a.parentSC->superClusterRef()->energy() :
1120  a.primaryGSFs[0].first->GsftrackRef()->pt() );
1121  const double b_en = ( b.parentSC ?
1122  b.parentSC->superClusterRef()->energy() :
1123  b.primaryGSFs[0].first->GsftrackRef()->pt() );
1124  return a_en < b_en;
1125  });
1126  */
1127  _refinableObjects.insert(_refinableObjects.end(),fromGSF);
1128  } // end loop on GSF elements of block
1129 }
1130 
1131  bool PFEGammaAlgo::
1132  unwrapSuperCluster(const PFSCElement* thesc,
1133  std::vector<PFClusterFlaggedElement>& ecalclusters,
1134  ClusterMap& ecal2ps) {
1135  ecalclusters.clear();
1136  ecal2ps.clear();
1137  LOGVERB("PFEGammaAlgo")
1138  << "Pointer to SC element: 0x"
1139  << std::hex << thesc << std::dec << std::endl
1140  << "cleared ecalclusters and ecal2ps!" << std::endl;
1141  auto ecalbegin = _splayedblock[reco::PFBlockElement::ECAL].begin();
1142  auto ecalend = _splayedblock[reco::PFBlockElement::ECAL].end();
1143  if( ecalbegin == ecalend ) {
1144  LOGERR("PFEGammaAlgo::unwrapSuperCluster()")
1145  << "There are no ECAL elements in a block with imported SC!"
1146  << " This is a bug we should fix this!"
1147  << std::endl;
1148  return false;
1149  }
1150  reco::SuperClusterRef scref = thesc->superClusterRef();
1151  const bool is_pf_sc = thesc->fromPFSuperCluster();
1152  if( !(scref.isAvailable() && scref.isNonnull()) ) {
1153  throw cms::Exception("PFEGammaAlgo::unwrapSuperCluster()")
1154  << "SuperCluster pointed to by block element is null!"
1155  << std::endl;
1156  }
1157  LOGDRESSED("PFEGammaAlgo")
1158  << "Got a valid super cluster ref! 0x"
1159  << std::hex << scref.get() << std::dec << std::endl;
1160  const size_t nscclusters = scref->clustersSize();
1161  const size_t nscpsclusters = scref->preshowerClustersSize();
1162  size_t npfpsclusters = 0;
1163  size_t npfclusters = 0;
1164  LOGDRESSED("PFEGammaAlgo")
1165  << "Precalculated cluster multiplicities: "
1166  << nscclusters << ' ' << nscpsclusters << std::endl;
1167  NotCloserToOther<reco::PFBlockElement::SC,reco::PFBlockElement::ECAL>
1168  ecalClustersInSC(_currentblock,_currentlinks,thesc);
1169  auto firstnotinsc = std::partition(ecalbegin,ecalend,ecalClustersInSC);
1170  //reset the begin and end iterators
1171  ecalbegin = _splayedblock[reco::PFBlockElement::ECAL].begin();
1172  ecalend = _splayedblock[reco::PFBlockElement::ECAL].end();
1173 
1174  //get list of associated clusters by det id and energy matching
1175  //(only needed when using non-pf supercluster)
1176  std::vector<const ClusterElement*> safePFClusters = is_pf_sc ? std::vector<const ClusterElement*>() : getSCAssociatedECALsSafe(scref,_splayedblock[reco::PFBlockElement::ECAL]);
1177 
1178  if( firstnotinsc == ecalbegin ) {
1179  LOGERR("PFEGammaAlgo::unwrapSuperCluster()")
1180  << "No associated block elements to SuperCluster!"
1181  << " This is a bug we should fix!"
1182  << std::endl;
1183  return false;
1184  }
1185  npfclusters = std::distance(ecalbegin,firstnotinsc);
1186  // ensure we have found the correct number of PF ecal clusters in the case
1187  // that this is a PF supercluster, otherwise all bets are off
1188  if( is_pf_sc && nscclusters != npfclusters ) {
1189  std::stringstream sc_err;
1190  thesc->Dump(sc_err,"\t");
1191  throw cms::Exception("PFEGammaAlgo::unwrapSuperCluster()")
1192  << "The number of found ecal elements ("
1193  << nscclusters << ") in block is not the same as"
1194  << " the number of ecal PF clusters reported by the PFSuperCluster"
1195  << " itself (" << npfclusters
1196  << ")! This should not happen!" << std::endl
1197  << sc_err.str() << std::endl;
1198  }
1199  for( auto ecalitr = ecalbegin; ecalitr != firstnotinsc; ++ecalitr ) {
1200  const PFClusterElement* elemascluster =
1201  docast(const PFClusterElement*,ecalitr->first);
1202 
1203  // reject clusters that really shouldn't be associated to the SC
1204  // (only needed when using non-pf-supercluster)
1205  if(!is_pf_sc && std::find(safePFClusters.begin(),safePFClusters.end(),elemascluster) ==
1206  safePFClusters.end() ) continue;
1207 
1208  //add cluster
1209  ecalclusters.push_back(std::make_pair(elemascluster,true));
1210  //mark cluster as used
1211  ecalitr->second = false;
1212 
1213  // process the ES elements
1214  // auto is a pair<Iterator,bool> here, bool is false when placing fails
1215  auto emplaceresult = ecal2ps.emplace(elemascluster,
1216  ClusterMap::mapped_type());
1217  if( !emplaceresult.second ) {
1218  std::stringstream clus_err;
1219  elemascluster->Dump(clus_err,"\t");
1220  throw cms::Exception("PFEGammaAlgo::unwrapSuperCluster()")
1221  << "List of pointers to ECAL block elements contains non-unique items!"
1222  << " This is very bad!" << std::endl
1223  << "cluster ptr = 0x" << std::hex << elemascluster << std::dec
1224  << std::endl << clus_err.str() << std::endl;
1225  }
1226  ClusterMap::mapped_type& eslist = emplaceresult.first->second;
1227  npfpsclusters += attachPSClusters(elemascluster,eslist);
1228  } // loop over ecal elements
1229 
1230  /*
1231  if( is_pf_sc && nscpsclusters != npfpsclusters) {
1232  std::stringstream sc_err;
1233  thesc->Dump(sc_err,"\t");
1234  throw cms::Exception("PFEGammaAlgo::unwrapSuperCluster()")
1235  << "The number of found PF preshower elements ("
1236  << npfpsclusters << ") in block is not the same as"
1237  << " the number of preshower clusters reported by the PFSuperCluster"
1238  << " itself (" << nscpsclusters << ")! This should not happen!"
1239  << std::endl
1240  << sc_err.str() << std::endl;
1241  }
1242  */
1243 
1244  LOGDRESSED("PFEGammaAlgo")
1245  << " Unwrapped SC has " << npfclusters << " ECAL sub-clusters"
1246  << " and " << npfpsclusters << " PreShower layers 1 & 2 clusters!"
1247  << std::endl;
1248  return true;
1249  }
1250 
1251 
1252 
1253  int PFEGammaAlgo::attachPSClusters(const ClusterElement* ecalclus,
1254  ClusterMap::mapped_type& eslist) {
1255  if( ecalclus->clusterRef()->layer() == PFLayer::ECAL_BARREL ) return 0;
1256  edm::Ptr<reco::PFCluster> clusptr = refToPtr(ecalclus->clusterRef());
1257  EEtoPSElement ecalkey = std::make_pair(clusptr.key(),clusptr);
1258  auto assc_ps = std::equal_range(eetops_->cbegin(),
1259  eetops_->cend(),
1260  ecalkey,
1261  comparePSMapByKey);
1262  for( const auto& ps1 : _splayedblock[reco::PFBlockElement::PS1] ) {
1263  edm::Ptr<reco::PFCluster> temp = refToPtr(ps1.first->clusterRef());
1264  for( auto pscl = assc_ps.first; pscl != assc_ps.second; ++pscl ) {
1265  if( pscl->second == temp ) {
1266  const ClusterElement* pstemp =
1267  docast(const ClusterElement*,ps1.first);
1268  eslist.push_back( PFClusterFlaggedElement(pstemp,true) );
1269  }
1270  }
1271  }
1272  for( const auto& ps2 : _splayedblock[reco::PFBlockElement::PS2] ) {
1273  edm::Ptr<reco::PFCluster> temp = refToPtr(ps2.first->clusterRef());
1274  for( auto pscl = assc_ps.first; pscl != assc_ps.second; ++pscl ) {
1275  if( pscl->second == temp ) {
1276  const ClusterElement* pstemp =
1277  docast(const ClusterElement*,ps2.first);
1278  eslist.push_back( PFClusterFlaggedElement(pstemp,true) );
1279  }
1280  }
1281  }
1282  return eslist.size();
1283  }
1284 
1285  void PFEGammaAlgo::dumpCurrentRefinableObjects() const {
1286  #ifdef PFLOW_DEBUG
1287  edm::LogVerbatim("PFEGammaAlgo")
1288  //<< "Dumping current block: " << std::endl << *_currentblock << std::endl
1289  << "Dumping " << _refinableObjects.size()
1290  << " refinable objects for this block: " << std::endl;
1291  for( const auto& ro : _refinableObjects ) {
1292  std::stringstream info;
1293  info << "Refinable Object:" << std::endl;
1294  if( ro.parentSC ) {
1295  info << "\tSuperCluster element attached to object:" << std::endl
1296  << '\t';
1297  ro.parentSC->Dump(info,"\t");
1298  info << std::endl;
1299  }
1300  if( ro.electronSeed.isNonnull() ) {
1301  info << "\tGSF element attached to object:" << std::endl;
1302  ro.primaryGSFs.front().first->Dump(info,"\t");
1303  info << std::endl;
1304  info << "firstBrem : " << ro.firstBrem
1305  << " lateBrem : " << ro.lateBrem
1306  << " nBrems with cluster : " << ro.nBremsWithClusters
1307  << std::endl;;
1308  if( ro.electronClusters.size() && ro.electronClusters[0] ) {
1309  info << "electron cluster : ";
1310  ro.electronClusters[0]->Dump(info,"\t");
1311  info << std::endl;
1312  } else {
1313  info << " no electron cluster." << std::endl;
1314  }
1315  }
1316  if( ro.primaryKFs.size() ) {
1317  info << "\tPrimary KF tracks attached to object: " << std::endl;
1318  for( const auto& kf : ro.primaryKFs ) {
1319  kf.first->Dump(info,"\t");
1320  info << std::endl;
1321  }
1322  }
1323  if( ro.secondaryKFs.size() ) {
1324  info << "\tSecondary KF tracks attached to object: " << std::endl;
1325  for( const auto& kf : ro.secondaryKFs ) {
1326  kf.first->Dump(info,"\t");
1327  info << std::endl;
1328  }
1329  }
1330  if( ro.brems.size() ) {
1331  info << "\tBrem tangents attached to object: " << std::endl;
1332  for( const auto& brem : ro.brems ) {
1333  brem.first->Dump(info,"\t");
1334  info << std::endl;
1335  }
1336  }
1337  if( ro.ecalclusters.size() ) {
1338  info << "\tECAL clusters attached to object: " << std::endl;
1339  for( const auto& clus : ro.ecalclusters ) {
1340  clus.first->Dump(info,"\t");
1341  info << std::endl;
1342  if( ro.ecal2ps.find(clus.first) != ro.ecal2ps.end() ) {
1343  for( const auto& psclus : ro.ecal2ps.at(clus.first) ) {
1344  info << "\t\t Attached PS Cluster: ";
1345  psclus.first->Dump(info,"");
1346  info << std::endl;
1347  }
1348  }
1349  }
1350  }
1351  edm::LogVerbatim("PFEGammaAlgo") << info.str();
1352  }
1353  #endif
1354  }
1355 
1356  // look through our KF tracks in this block and match
1357  void PFEGammaAlgo::
1358  removeOrLinkECALClustersToKFTracks() {
1359  if( !_splayedblock[reco::PFBlockElement::ECAL].size() ||
1360  !_splayedblock[reco::PFBlockElement::TRACK].size() ) return;
1361  std::multimap<double, unsigned> matchedGSFs, matchedECALs;
1362  for( auto& kftrack : _splayedblock[reco::PFBlockElement::TRACK] ) {
1363  matchedGSFs.clear();
1364  _currentblock->associatedElements(kftrack.first->index(), _currentlinks,
1365  matchedGSFs,
1366  reco::PFBlockElement::GSF,
1367  reco::PFBlock::LINKTEST_ALL);
1368  if( !matchedGSFs.size() ) { // only run this is we aren't associated to GSF
1369  LesserByDistance closestTrackToECAL(_currentblock,_currentlinks,
1370  &kftrack);
1371  auto ecalbegin = _splayedblock[reco::PFBlockElement::ECAL].begin();
1372  auto ecalend = _splayedblock[reco::PFBlockElement::ECAL].end();
1373  std::partial_sort(ecalbegin,ecalbegin+1,ecalend,closestTrackToECAL);
1374  PFFlaggedElement& closestECAL =
1375  _splayedblock[reco::PFBlockElement::ECAL].front();
1376  const float dist = _currentblock->dist(kftrack.first->index(),
1377  closestECAL.first->index(),
1378  _currentlinks,
1379  reco::PFBlock::LINKTEST_ALL);
1380  bool inSC = false;
1381  for( auto& sc : _splayedblock[reco::PFBlockElement::SC] ) {
1382  float dist_sc = _currentblock->dist(sc.first->index(),
1383  closestECAL.first->index(),
1384  _currentlinks,
1385  reco::PFBlock::LINKTEST_ALL);
1386  if( dist_sc != -1.0f) { inSC = true; break; }
1387  }
1388 
1389  if( dist != -1.0f && closestECAL.second ) {
1390  bool gsflinked = false;
1391  // check that this cluster is not associated to a GSF track
1392  for(const auto& gsfflag : _splayedblock[reco::PFBlockElement::GSF]) {
1393  const reco::PFBlockElementGsfTrack* elemasgsf =
1394  docast(const reco::PFBlockElementGsfTrack*,gsfflag.first);
1395  if(elemasgsf->trackType(reco::PFBlockElement::T_FROM_GAMMACONV)) {
1396  continue; // keep clusters that have a found conversion GSF near
1397  }
1398  matchedECALs.clear();
1399  _currentblock->associatedElements(elemasgsf->index(), _currentlinks,
1400  matchedECALs,
1401  reco::PFBlockElement::ECAL,
1402  reco::PFBlock::LINKTEST_ALL);
1403  if( matchedECALs.size() ) {
1404  if( matchedECALs.begin()->second == closestECAL.first->index() ) {
1405  gsflinked = true;
1406  break;
1407  }
1408  }
1409  } // loop over primary GSF tracks
1410  if( !gsflinked && !inSC) {
1411  // determine if we should remove the matched cluster
1412  const reco::PFBlockElementTrack * kfEle =
1413  docast(const reco::PFBlockElementTrack*,kftrack.first);
1414  const reco::TrackRef trackref = kfEle->trackRef();
1415  const unsigned Algo = trackref->algo();
1416  const int nexhits =
1417  trackref->trackerExpectedHitsInner().numberOfLostHits();
1418  bool fromprimaryvertex = false;
1419  for( auto vtxtks = cfg_.primaryVtx->tracks_begin();
1420  vtxtks != cfg_.primaryVtx->tracks_end(); ++ vtxtks ) {
1421  if( trackref == vtxtks->castTo<reco::TrackRef>() ) {
1422  fromprimaryvertex = true;
1423  break;
1424  }
1425  }// loop over tracks in primary vertex
1426  // if associated to good non-GSF matched track remove this cluster
1427  if( Algo < 9 && nexhits == 0 && fromprimaryvertex ) {
1428  closestECAL.second = false;
1429  } else { // otherwise associate the cluster and KF track
1430  _recoveredlinks.push_back( ElementMap::value_type(closestECAL.first,kftrack.first) );
1431  _recoveredlinks.push_back( ElementMap::value_type(kftrack.first,closestECAL.first) );
1432  }
1433  }
1434  } // found a good closest ECAL match
1435  } // no GSF track matched to KF
1436  } // loop over KF elements
1437  }
1438 
1439  void PFEGammaAlgo::
1440  mergeROsByAnyLink(std::list<PFEGammaAlgo::ProtoEGObject>& ROs) {
1441  if( ROs.size() < 2 ) return; // nothing to do with one or zero ROs
1442  bool check_for_merge = true;
1443  while( check_for_merge ) {
1444  // bugfix for early termination merging loop (15 April 2014)
1445  // check all pairwise combinations in the list
1446  // if one has a merge shuffle it to the front of the list
1447  // if there are no merges left to do we can terminate
1448  for( auto it1 = ROs.begin(); it1 != ROs.end(); ++it1 ) {
1449  TestIfROMergableByLink mergeTest(*it1);
1450  auto find_start = it1; ++find_start;
1451  auto has_merge = std::find_if(find_start,ROs.end(),mergeTest);
1452  if( has_merge != ROs.end() && it1 != ROs.begin() ) {
1453  std::swap(*(ROs.begin()),*it1);
1454  break;
1455  }
1456  }// ensure mergables are shuffled to the front
1457  ProtoEGObject& thefront = ROs.front();
1458  TestIfROMergableByLink mergeTest(thefront);
1459  auto mergestart = ROs.begin(); ++mergestart;
1460  auto nomerge = std::partition(mergestart,ROs.end(),mergeTest);
1461  if( nomerge != mergestart ) {
1462  LOGDRESSED("PFEGammaAlgo::mergeROsByAnyLink()")
1463  << "Found objects " << std::distance(mergestart,nomerge)
1464  << " to merge by links to the front!" << std::endl;
1465  for( auto roToMerge = mergestart; roToMerge != nomerge; ++roToMerge) {
1466  thefront.ecalclusters.insert(thefront.ecalclusters.end(),
1467  roToMerge->ecalclusters.begin(),
1468  roToMerge->ecalclusters.end());
1469  thefront.ecal2ps.insert(roToMerge->ecal2ps.begin(),
1470  roToMerge->ecal2ps.end());
1471  thefront.secondaryKFs.insert(thefront.secondaryKFs.end(),
1472  roToMerge->secondaryKFs.begin(),
1473  roToMerge->secondaryKFs.end());
1474 
1475  thefront.localMap.insert(thefront.localMap.end(),
1476  roToMerge->localMap.begin(),
1477  roToMerge->localMap.end());
1478  // TO FIX -> use best (E_gsf - E_clustersum)/E_GSF
1479  if( !thefront.parentSC && roToMerge->parentSC ) {
1480  thefront.parentSC = roToMerge->parentSC;
1481  }
1482  if( thefront.electronSeed.isNull() &&
1483  roToMerge->electronSeed.isNonnull() ) {
1484  thefront.electronSeed = roToMerge->electronSeed;
1485  thefront.primaryGSFs.insert(thefront.primaryGSFs.end(),
1486  roToMerge->primaryGSFs.begin(),
1487  roToMerge->primaryGSFs.end());
1488  thefront.primaryKFs.insert(thefront.primaryKFs.end(),
1489  roToMerge->primaryKFs.begin(),
1490  roToMerge->primaryKFs.end());
1491  thefront.brems.insert(thefront.brems.end(),
1492  roToMerge->brems.begin(),
1493  roToMerge->brems.end());
1494  thefront.electronClusters = roToMerge->electronClusters;
1495  thefront.nBremsWithClusters = roToMerge->nBremsWithClusters;
1496  thefront.firstBrem = roToMerge->firstBrem;
1497  thefront.lateBrem = roToMerge->lateBrem;
1498  } else if ( thefront.electronSeed.isNonnull() &&
1499  roToMerge->electronSeed.isNonnull()) {
1500  LOGWARN("PFEGammaAlgo::mergeROsByAnyLink")
1501  << "Need to implement proper merging of two gsf candidates!"
1502  << std::endl;
1503  }
1504  }
1505  ROs.erase(mergestart,nomerge);
1506  // put the merged element in the back of the cleaned list
1507  ROs.push_back(ROs.front());
1508  ROs.pop_front();
1509  } else {
1510  check_for_merge = false;
1511  }
1512  }
1513  LOGDRESSED("PFEGammaAlgo::mergeROsByAnyLink()")
1514  << "After merging by links there are: " << ROs.size()
1515  << " refinable EGamma objects!" << std::endl;
1516  }
1517 
1518 // pull in KF tracks associated to the RO but not closer to another
1519 // NB: in initializeProtoCands() we forced the GSF tracks not to be
1520 // from a conversion, but we will leave a protection here just in
1521 // case things change in the future
1522 void PFEGammaAlgo::
1523 linkRefinableObjectGSFTracksToKFs(ProtoEGObject& RO) {
1524  constexpr reco::PFBlockElement::TrackType convType =
1525  reco::PFBlockElement::T_FROM_GAMMACONV;
1526  if( !_splayedblock[reco::PFBlockElement::TRACK].size() ) return;
1527  auto KFbegin = _splayedblock[reco::PFBlockElement::TRACK].begin();
1528  auto KFend = _splayedblock[reco::PFBlockElement::TRACK].end();
1529  for( auto& gsfflagged : RO.primaryGSFs ) {
1530  const PFGSFElement* seedtk = gsfflagged.first;
1531  // don't process SC-only ROs or secondary seeded ROs
1532  if( RO.electronSeed.isNull() || seedtk->trackType(convType) ) continue;
1533  NotCloserToOther<reco::PFBlockElement::GSF,reco::PFBlockElement::TRACK>
1534  gsfTrackToKFs(_currentblock,_currentlinks,seedtk);
1535  // get KF tracks not closer to another and not already used
1536  auto notlinked = std::partition(KFbegin,KFend,gsfTrackToKFs);
1537  // attach tracks and set as used
1538  for( auto kft = KFbegin; kft != notlinked; ++kft ) {
1539  const PFKFElement* elemaskf =
1540  docast(const PFKFElement*,kft->first);
1541  // don't care about things that aren't primaries or directly
1542  // associated secondary tracks
1543  if( isPrimaryTrack(*elemaskf,*seedtk) &&
1544  !elemaskf->trackType(convType) ) {
1545  kft->second = false;
1546  RO.primaryKFs.push_back(std::make_pair(elemaskf,true));
1547  RO.localMap.push_back( ElementMap::value_type(seedtk,elemaskf) );
1548  RO.localMap.push_back( ElementMap::value_type(elemaskf,seedtk) );
1549  } else if ( elemaskf->trackType(convType) ) {
1550  kft->second = false;
1551  RO.secondaryKFs.push_back(std::make_pair(elemaskf,true));
1552  RO.localMap.push_back( ElementMap::value_type(seedtk,elemaskf) );
1553  RO.localMap.push_back( ElementMap::value_type(elemaskf,seedtk) );
1554  }
1555  }// loop on closest KFs not closer to other GSFs
1556  } // loop on GSF primaries on RO
1557 }
1558 
1559 void PFEGammaAlgo::
1560 linkRefinableObjectPrimaryKFsToSecondaryKFs(ProtoEGObject& RO) {
1561  constexpr reco::PFBlockElement::TrackType convType =
1562  reco::PFBlockElement::T_FROM_GAMMACONV;
1563  if( !_splayedblock[reco::PFBlockElement::TRACK].size() ) return;
1564  auto KFbegin = _splayedblock[reco::PFBlockElement::TRACK].begin();
1565  auto KFend = _splayedblock[reco::PFBlockElement::TRACK].end();
1566  for( auto& kfflagged : RO.primaryKFs ) {
1567  const PFKFElement* primkf = kfflagged.first;
1568  // don't process SC-only ROs or secondary seeded ROs
1569  if( primkf->trackType(convType) ) {
1570  throw cms::Exception("PFEGammaAlgo::linkRefinableObjectPrimaryKFsToSecondaryKFs()")
1571  << "A KF track from conversion has been assigned as a primary!!"
1572  << std::endl;
1573  }
1574  NotCloserToOther<reco::PFBlockElement::TRACK,reco::PFBlockElement::TRACK,true>
1575  kfTrackToKFs(_currentblock,_currentlinks,primkf);
1576  // get KF tracks not closer to another and not already used
1577  auto notlinked = std::partition(KFbegin,KFend,kfTrackToKFs);
1578  // attach tracks and set as used
1579  for( auto kft = KFbegin; kft != notlinked; ++kft ) {
1580  const PFKFElement* elemaskf =
1581  docast(const PFKFElement*,kft->first);
1582  // don't care about things that aren't primaries or directly
1583  // associated secondary tracks
1584  if( elemaskf->trackType(convType) ) {
1585  kft->second = false;
1586  RO.secondaryKFs.push_back(std::make_pair(elemaskf,true));
1587  RO.localMap.push_back( ElementMap::value_type(primkf,elemaskf) );
1588  RO.localMap.push_back( ElementMap::value_type(elemaskf,primkf) );
1589  }
1590  }// loop on closest KFs not closer to other KFs
1591  } // loop on KF primaries on RO
1592 }
1593 
1594 // try to associate the tracks to cluster elements which are not used
1595 void PFEGammaAlgo::
1596 linkRefinableObjectPrimaryGSFTrackToECAL(ProtoEGObject& RO) {
1597  if( !_splayedblock[reco::PFBlockElement::ECAL].size() ) {
1598  RO.electronClusters.push_back(NULL);
1599  return;
1600  }
1601  auto ECALbegin = _splayedblock[reco::PFBlockElement::ECAL].begin();
1602  auto ECALend = _splayedblock[reco::PFBlockElement::ECAL].end();
1603  for( auto& primgsf : RO.primaryGSFs ) {
1604  NotCloserToOther<reco::PFBlockElement::GSF,reco::PFBlockElement::ECAL>
1605  gsfTracksToECALs(_currentblock,_currentlinks,primgsf.first);
1606  CompatibleEoPOut eoverp_test(primgsf.first);
1607  // get set of matching ecals not already in SC
1608  auto notmatched_blk = std::partition(ECALbegin,ECALend,gsfTracksToECALs);
1609  notmatched_blk = std::partition(ECALbegin,notmatched_blk,eoverp_test);
1610  // get set of matching ecals already in the RO
1611  auto notmatched_sc = std::partition(RO.ecalclusters.begin(),
1612  RO.ecalclusters.end(),
1613  gsfTracksToECALs);
1614  notmatched_sc = std::partition(RO.ecalclusters.begin(),
1615  notmatched_sc,
1616  eoverp_test);
1617  // look inside the SC for the ECAL cluster
1618  for( auto ecal = RO.ecalclusters.begin(); ecal != notmatched_sc; ++ecal ) {
1619  const PFClusterElement* elemascluster =
1620  docast(const PFClusterElement*,ecal->first);
1621  PFClusterFlaggedElement temp(elemascluster,true);
1622  LOGDRESSED("PFEGammaAlgo::linkGSFTracktoECAL()")
1623  << "Found a cluster already in RO by GSF extrapolation"
1624  << " at ECAL surface!" << std::endl
1625  << *elemascluster << std::endl;
1626 
1627  RO.localMap.push_back(ElementMap::value_type(primgsf.first,temp.first));
1628  RO.localMap.push_back(ElementMap::value_type(temp.first,primgsf.first));
1629  }
1630  // look outside the SC for the ecal cluster
1631  for( auto ecal = ECALbegin; ecal != notmatched_blk; ++ecal ) {
1632  const PFClusterElement* elemascluster =
1633  docast(const PFClusterElement*,ecal->first);
1634  LOGDRESSED("PFEGammaAlgo::linkGSFTracktoECAL()")
1635  << "Found a cluster not already in RO by GSF extrapolation"
1636  << " at ECAL surface!" << std::endl
1637  << *elemascluster << std::endl;
1638  if( addPFClusterToROSafe(elemascluster,RO) ) {
1639  attachPSClusters(elemascluster,RO.ecal2ps[elemascluster]);
1640  RO.localMap.push_back(ElementMap::value_type(primgsf.first,elemascluster));
1641  RO.localMap.push_back(ElementMap::value_type(elemascluster,primgsf.first));
1642  ecal->second = false;
1643  }
1644  }
1645  }
1646 }
1647 
1648 // try to associate the tracks to cluster elements which are not used
1649 void PFEGammaAlgo::
1650 linkRefinableObjectPrimaryGSFTrackToHCAL(ProtoEGObject& RO) {
1651  if( !_splayedblock[reco::PFBlockElement::HCAL].size() ) return;
1652  auto HCALbegin = _splayedblock[reco::PFBlockElement::HCAL].begin();
1653  auto HCALend = _splayedblock[reco::PFBlockElement::HCAL].end();
1654  for( auto& primgsf : RO.primaryGSFs ) {
1655  NotCloserToOther<reco::PFBlockElement::GSF,reco::PFBlockElement::HCAL>
1656  gsfTracksToHCALs(_currentblock,_currentlinks,primgsf.first);
1657  CompatibleEoPOut eoverp_test(primgsf.first);
1658  auto notmatched = std::partition(HCALbegin,HCALend,gsfTracksToHCALs);
1659  for( auto hcal = HCALbegin; hcal != notmatched; ++hcal ) {
1660  const PFClusterElement* elemascluster =
1661  docast(const PFClusterElement*,hcal->first);
1662  PFClusterFlaggedElement temp(elemascluster,true);
1663  LOGDRESSED("PFEGammaAlgo::linkGSFTracktoECAL()")
1664  << "Found an HCAL cluster associated to GSF extrapolation"
1665  << std::endl;
1666  RO.hcalClusters.push_back(temp);
1667  RO.localMap.push_back( ElementMap::value_type(primgsf.first,temp.first) );
1668  RO.localMap.push_back( ElementMap::value_type(temp.first,primgsf.first) );
1669  hcal->second = false;
1670  }
1671  }
1672 }
1673 
1674 // try to associate the tracks to cluster elements which are not used
1675 void PFEGammaAlgo::
1676 linkRefinableObjectKFTracksToECAL(ProtoEGObject& RO) {
1677  if( !_splayedblock[reco::PFBlockElement::ECAL].size() ) return;
1678  for( auto& primkf : RO.primaryKFs ) linkKFTrackToECAL(primkf,RO);
1679  for( auto& secdkf : RO.secondaryKFs ) linkKFTrackToECAL(secdkf,RO);
1680 }
1681 
1682 void
1683 PFEGammaAlgo::linkKFTrackToECAL(const KFFlaggedElement& kfflagged,
1684  ProtoEGObject& RO) {
1685  std::vector<PFClusterFlaggedElement>& currentECAL = RO.ecalclusters;
1686  auto ECALbegin = _splayedblock[reco::PFBlockElement::ECAL].begin();
1687  auto ECALend = _splayedblock[reco::PFBlockElement::ECAL].end();
1688  NotCloserToOther<reco::PFBlockElement::TRACK,reco::PFBlockElement::ECAL>
1689  kfTrackToECALs(_currentblock,_currentlinks,kfflagged.first);
1690  NotCloserToOther<reco::PFBlockElement::GSF,reco::PFBlockElement::ECAL>
1691  kfTrackGSFToECALs(_currentblock,_currentlinks,kfflagged.first);
1692  //get the ECAL elements not used and not closer to another KF
1693  auto notmatched_sc = std::partition(currentECAL.begin(),
1694  currentECAL.end(),
1695  kfTrackToECALs);
1696  //get subset ECAL elements not used or closer to another GSF of any type
1697  notmatched_sc = std::partition(currentECAL.begin(),
1698  notmatched_sc,
1699  kfTrackGSFToECALs);
1700  for( auto ecalitr = currentECAL.begin(); ecalitr != notmatched_sc;
1701  ++ecalitr ) {
1702  const PFClusterElement* elemascluster =
1703  docast(const PFClusterElement*,ecalitr->first);
1704  PFClusterFlaggedElement flaggedclus(elemascluster,true);
1705 
1706  LOGDRESSED("PFEGammaAlgo::linkKFTracktoECAL()")
1707  << "Found a cluster already in RO by KF extrapolation"
1708  << " at ECAL surface!" << std::endl
1709  << *elemascluster << std::endl;
1710  RO.localMap.push_back(ElementMap::value_type(elemascluster,
1711  kfflagged.first));
1712  RO.localMap.push_back(ElementMap::value_type(kfflagged.first,
1713  elemascluster));
1714  }
1715  //get the ECAL elements not used and not closer to another KF
1716  auto notmatched_blk = std::partition(ECALbegin,ECALend,kfTrackToECALs);
1717  //get subset ECAL elements not used or closer to another GSF of any type
1718  notmatched_blk = std::partition(ECALbegin,notmatched_blk,kfTrackGSFToECALs);
1719  for( auto ecalitr = ECALbegin; ecalitr != notmatched_blk; ++ecalitr ) {
1720  const PFClusterElement* elemascluster =
1721  docast(const PFClusterElement*,ecalitr->first);
1722  if( addPFClusterToROSafe(elemascluster,RO) ) {
1723  attachPSClusters(elemascluster,RO.ecal2ps[elemascluster]);
1724  ecalitr->second = false;
1725 
1726  LOGDRESSED("PFEGammaAlgo::linkKFTracktoECAL()")
1727  << "Found a cluster not in RO by KF extrapolation"
1728  << " at ECAL surface!" << std::endl
1729  << *elemascluster << std::endl;
1730  RO.localMap.push_back(ElementMap::value_type(elemascluster,
1731  kfflagged.first));
1732  RO.localMap.push_back( ElementMap::value_type(kfflagged.first,
1733  elemascluster));
1734  }
1735  }
1736 }
1737 
1738 void PFEGammaAlgo::
1739 linkRefinableObjectBremTangentsToECAL(ProtoEGObject& RO) {
1740  if( !RO.brems.size() ) return;
1741  int FirstBrem = -1;
1742  int TrajPos = -1;
1743  int lastBremTrajPos = -1;
1744  for( auto& bremflagged : RO.brems ) {
1745  bool has_clusters = false;
1746  TrajPos = (bremflagged.first->indTrajPoint())-2;
1747  auto ECALbegin = _splayedblock[reco::PFBlockElement::ECAL].begin();
1748  auto ECALend = _splayedblock[reco::PFBlockElement::ECAL].end();
1749  NotCloserToOther<reco::PFBlockElement::BREM,reco::PFBlockElement::ECAL>
1750  BremToECALs(_currentblock,_currentlinks,bremflagged.first);
1751  // check for late brem using clusters already in the SC
1752  auto RSCBegin = RO.ecalclusters.begin();
1753  auto RSCEnd = RO.ecalclusters.end();
1754  auto notmatched_rsc = std::partition(RSCBegin,RSCEnd,BremToECALs);
1755  for( auto ecal = RSCBegin; ecal != notmatched_rsc; ++ecal ) {
1756  float deta =
1757  std::abs( ecal->first->clusterRef()->positionREP().eta() -
1758  bremflagged.first->positionAtECALEntrance().eta() );
1759  if( deta < 0.015 ) {
1760  has_clusters = true;
1761  if( lastBremTrajPos == -1 || lastBremTrajPos < TrajPos ) {
1762  lastBremTrajPos = TrajPos;
1763  }
1764  if( FirstBrem == -1 || TrajPos < FirstBrem ) { // set brem information
1765  FirstBrem = TrajPos;
1766  RO.firstBrem = TrajPos;
1767  }
1768  LOGDRESSED("PFEGammaAlgo::linkBremToECAL()")
1769  << "Found a cluster already in SC linked to brem extrapolation"
1770  << " at ECAL surface!" << std::endl;
1771  RO.localMap.push_back( ElementMap::value_type(ecal->first,bremflagged.first) );
1772  RO.localMap.push_back( ElementMap::value_type(bremflagged.first,ecal->first) );
1773  }
1774  }
1775  // grab new clusters from the block (ensured to not be late brem)
1776  auto notmatched_block = std::partition(ECALbegin,ECALend,BremToECALs);
1777  for( auto ecal = ECALbegin; ecal != notmatched_block; ++ecal ) {
1778  float deta =
1779  std::abs( ecal->first->clusterRef()->positionREP().eta() -
1780  bremflagged.first->positionAtECALEntrance().eta() );
1781  if( deta < 0.015 ) {
1782  has_clusters = true;
1783  if( lastBremTrajPos == -1 || lastBremTrajPos < TrajPos ) {
1784  lastBremTrajPos = TrajPos;
1785  }
1786  if( FirstBrem == -1 || TrajPos < FirstBrem ) { // set brem information
1787 
1788  FirstBrem = TrajPos;
1789  RO.firstBrem = TrajPos;
1790  }
1791  const PFClusterElement* elemasclus =
1792  docast(const PFClusterElement*,ecal->first);
1793  if( addPFClusterToROSafe(elemasclus,RO) ) {
1794  attachPSClusters(elemasclus,RO.ecal2ps[elemasclus]);
1795 
1796  RO.localMap.push_back( ElementMap::value_type(ecal->first,bremflagged.first) );
1797  RO.localMap.push_back( ElementMap::value_type(bremflagged.first,ecal->first) );
1798  ecal->second = false;
1799  LOGDRESSED("PFEGammaAlgo::linkBremToECAL()")
1800  << "Found a cluster not already associated by brem extrapolation"
1801  << " at ECAL surface!" << std::endl;
1802  }
1803 
1804  }
1805  }
1806  if(has_clusters) {
1807  if( RO.nBremsWithClusters == -1 ) RO.nBremsWithClusters = 0;
1808  ++RO.nBremsWithClusters;
1809  }
1810  }
1811 }
1812 
1813 void PFEGammaAlgo::
1814 linkRefinableObjectConvSecondaryKFsToSecondaryKFs(ProtoEGObject& RO) {
1815  IsConversionTrack<reco::PFBlockElementTrack> isConvKf;
1816  auto KFbegin = _splayedblock[reco::PFBlockElement::TRACK].begin();
1817  auto KFend = _splayedblock[reco::PFBlockElement::TRACK].end();
1818  auto BeginROskfs = RO.secondaryKFs.begin();
1819  auto EndROskfs = RO.secondaryKFs.end();
1820  auto ronotconv = std::partition(BeginROskfs,EndROskfs,isConvKf);
1821  size_t convkfs_end = std::distance(BeginROskfs,ronotconv);
1822  for( size_t idx = 0; idx < convkfs_end; ++idx ) {
1823  const PFKFFlaggedElement ro_skf = RO.secondaryKFs[idx];
1824  NotCloserToOther<reco::PFBlockElement::TRACK,
1825  reco::PFBlockElement::TRACK,
1826  true>
1827  TracksToTracks(_currentblock,_currentlinks, ro_skf.first);
1828  auto notmatched = std::partition(KFbegin,KFend,TracksToTracks);
1829  notmatched = std::partition(KFbegin,notmatched,isConvKf);
1830  for( auto kf = KFbegin; kf != notmatched; ++kf ) {
1831  const reco::PFBlockElementTrack* elemaskf =
1832  docast(const reco::PFBlockElementTrack*,kf->first);
1833  RO.secondaryKFs.push_back( std::make_pair(elemaskf,true) );
1834  RO.localMap.push_back( ElementMap::value_type(ro_skf.first,kf->first) );
1835  RO.localMap.push_back( ElementMap::value_type(kf->first,ro_skf.first) );
1836  kf->second = false;
1837  }
1838  }
1839 }
1840 
1841 void PFEGammaAlgo::
1842 linkRefinableObjectECALToSingleLegConv(ProtoEGObject& RO) {
1843  IsConversionTrack<reco::PFBlockElementTrack> isConvKf;
1844  auto KFbegin = _splayedblock[reco::PFBlockElement::TRACK].begin();
1845  auto KFend = _splayedblock[reco::PFBlockElement::TRACK].end();
1846  for( auto& ecal : RO.ecalclusters ) {
1847  NotCloserToOther<reco::PFBlockElement::ECAL,
1848  reco::PFBlockElement::TRACK,
1849  true>
1850  ECALToTracks(_currentblock,_currentlinks,ecal.first);
1851  auto notmatchedkf = std::partition(KFbegin,KFend,ECALToTracks);
1852  auto notconvkf = std::partition(KFbegin,notmatchedkf,isConvKf);
1853  // add identified KF conversion tracks
1854  for( auto kf = KFbegin; kf != notconvkf; ++kf ) {
1855  const reco::PFBlockElementTrack* elemaskf =
1856  docast(const reco::PFBlockElementTrack*,kf->first);
1857  RO.secondaryKFs.push_back( std::make_pair(elemaskf,true) );
1858  RO.localMap.push_back( ElementMap::value_type(ecal.first,elemaskf) );
1859  RO.localMap.push_back( ElementMap::value_type(elemaskf,ecal.first) );
1860  kf->second = false;
1861  }
1862  // go through non-conv-identified kfs and check MVA to add conversions
1863  for( auto kf = notconvkf; kf != notmatchedkf; ++kf ) {
1864  float mvaval = EvaluateSingleLegMVA(_currentblock, *cfg_.primaryVtx,
1865  kf->first->index());
1866  if(mvaval > cfg_.mvaConvCut) {
1867  const reco::PFBlockElementTrack* elemaskf =
1868  docast(const reco::PFBlockElementTrack*,kf->first);
1869  RO.secondaryKFs.push_back( std::make_pair(elemaskf,true) );
1870  RO.localMap.push_back( ElementMap::value_type(ecal.first,elemaskf) );
1871  RO.localMap.push_back( ElementMap::value_type(elemaskf,ecal.first) );
1872  kf->second = false;
1873 
1874  RO.singleLegConversionMvaMap.insert(std::make_pair(elemaskf, mvaval));
1875  }
1876  }
1877  }
1878 }
1879 
1880 void PFEGammaAlgo::
1881 linkRefinableObjectSecondaryKFsToECAL(ProtoEGObject& RO) {
1882  auto ECALbegin = _splayedblock[reco::PFBlockElement::ECAL].begin();
1883  auto ECALend = _splayedblock[reco::PFBlockElement::ECAL].end();
1884  for( auto& skf : RO.secondaryKFs ) {
1885  NotCloserToOther<reco::PFBlockElement::TRACK,
1886  reco::PFBlockElement::ECAL,
1887  false>
1888  TracksToECALwithCut(_currentblock,_currentlinks,skf.first,1.5f);
1889  auto notmatched = std::partition(ECALbegin,ECALend,TracksToECALwithCut);
1890  for( auto ecal = ECALbegin; ecal != notmatched; ++ecal ) {
1891  const reco::PFBlockElementCluster* elemascluster =
1892  docast(const reco::PFBlockElementCluster*,ecal->first);
1893  if( addPFClusterToROSafe(elemascluster,RO) ) {
1894  attachPSClusters(elemascluster,RO.ecal2ps[elemascluster]);
1895  RO.localMap.push_back(ElementMap::value_type(skf.first,elemascluster));
1896  RO.localMap.push_back(ElementMap::value_type(elemascluster,skf.first));
1897  ecal->second = false;
1898  }
1899  }
1900  }
1901 }
1902 
1903 void PFEGammaAlgo::
1904 fillPFCandidates(const std::list<PFEGammaAlgo::ProtoEGObject>& ROs,
1905  reco::PFCandidateCollection& egcands,
1906  reco::PFCandidateEGammaExtraCollection& egxs) {
1907  // reset output collections
1908  egcands.clear();
1909  egxs.clear();
1910  refinedscs_.clear();
1911  egcands.reserve(ROs.size());
1912  egxs.reserve(ROs.size());
1913  refinedscs_.reserve(ROs.size());
1914  for( auto& RO : ROs ) {
1915  if( RO.ecalclusters.size() == 0 &&
1916  !cfg_.produceEGCandsWithNoSuperCluster ) continue;
1917 
1918  reco::PFCandidate cand;
1919  reco::PFCandidateEGammaExtra xtra;
1920  if( RO.primaryGSFs.size() || RO.primaryKFs.size() ) {
1921  cand.setPdgId(-11); // anything with a primary track is an electron
1922  } else {
1923  cand.setPdgId(22); // anything with no primary track is a photon
1924  }
1925  if( RO.primaryKFs.size() ) {
1926  cand.setCharge(RO.primaryKFs[0].first->trackRef()->charge());
1927  xtra.setKfTrackRef(RO.primaryKFs[0].first->trackRef());
1928  cand.setTrackRef(RO.primaryKFs[0].first->trackRef());
1929  cand.addElementInBlock(_currentblock,RO.primaryKFs[0].first->index());
1930  }
1931  if( RO.primaryGSFs.size() ) {
1932  cand.setCharge(RO.primaryGSFs[0].first->GsftrackRef()->chargeMode());
1933  xtra.setGsfTrackRef(RO.primaryGSFs[0].first->GsftrackRef());
1934  cand.setGsfTrackRef(RO.primaryGSFs[0].first->GsftrackRef());
1935  cand.addElementInBlock(_currentblock,RO.primaryGSFs[0].first->index());
1936  }
1937  if( RO.parentSC ) {
1938  xtra.setSuperClusterPFECALRef(RO.parentSC->superClusterRef());
1939  // we'll set to the refined supercluster back up in the producer
1940  cand.setSuperClusterRef(RO.parentSC->superClusterRef());
1941  xtra.setSuperClusterRef(RO.parentSC->superClusterRef());
1942  cand.addElementInBlock(_currentblock,RO.parentSC->index());
1943  }
1944  // add brems
1945  for( const auto& bremflagged : RO.brems ) {
1946  const PFBremElement* brem = bremflagged.first;
1947  cand.addElementInBlock(_currentblock,brem->index());
1948  }
1949  // add clusters and ps clusters
1950  for( const auto& ecal : RO.ecalclusters ) {
1951  const PFClusterElement* clus = ecal.first;
1952  cand.addElementInBlock(_currentblock,clus->index());
1953  for( auto& ps : RO.ecal2ps.at(clus) ) {
1954  const PFClusterElement* psclus = ps.first;
1955  cand.addElementInBlock(_currentblock,psclus->index());
1956  }
1957  }
1958  // add secondary tracks
1959  for( const auto& secdkf : RO.secondaryKFs ) {
1960  const PFKFElement* kf = secdkf.first;
1961  cand.addElementInBlock(_currentblock,kf->index());
1962  const reco::ConversionRefVector& convrefs = kf->convRefs();
1963  bool no_conv_ref = true;
1964  for( const auto& convref : convrefs ) {
1965  if( convref.isNonnull() && convref.isAvailable() ) {
1966  xtra.addConversionRef(convref);
1967  no_conv_ref = false;
1968  }
1969  }
1970  if( no_conv_ref ) {
1971  //single leg conversions
1972 
1973  //look for stored mva value in map or else recompute
1974  const auto &mvavalmapped = RO.singleLegConversionMvaMap.find(kf);
1975  //FIXME: Abuse single mva value to store both provenance and single leg mva score
1976  //by storing 3.0 + mvaval
1977  float mvaval = mvavalmapped!=RO.singleLegConversionMvaMap.end() ? mvavalmapped->second : 3.0 + EvaluateSingleLegMVA(_currentblock, *cfg_.primaryVtx,
1978  kf->index());
1979 
1980  xtra.addSingleLegConvTrackRefMva(std::make_pair(kf->trackRef(),mvaval));
1981  }
1982  }
1983 
1984  // build the refined supercluster from those clusters left in the cand
1985  refinedscs_.push_back(buildRefinedSuperCluster(RO));
1986 
1987  const reco::SuperCluster& the_sc = refinedscs_.back();
1988  // with the refined SC in hand we build a naive candidate p4
1989  // and set the candidate ECAL position to either the barycenter of the
1990  // supercluster (if super-cluster present) or the seed of the
1991  // new SC generated by the EGAlgo
1992  const double scE = the_sc.energy();
1993  if( scE != 0.0 ) {
1994  const math::XYZPoint& seedPos = the_sc.seed()->position();
1995  math::XYZVector egDir = the_sc.position()-cfg_.primaryVtx->position();
1996  egDir = egDir.Unit();
1997  cand.setP4(math::XYZTLorentzVector(scE*egDir.x(),
1998  scE*egDir.y(),
1999  scE*egDir.z(),
2000  scE ));
2001  math::XYZPointF ecalPOS_f(seedPos.x(),seedPos.y(),seedPos.z());
2002  cand.setPositionAtECALEntrance(ecalPOS_f);
2003  cand.setEcalEnergy(the_sc.rawEnergy(),the_sc.energy());
2004  } else if ( cfg_.produceEGCandsWithNoSuperCluster &&
2005  RO.primaryGSFs.size() ) {
2006  const PFGSFElement* gsf = RO.primaryGSFs[0].first;
2007  reco::GsfTrackRef gref = gsf->GsftrackRef();
2008  math::XYZTLorentzVector p4(gref->pxMode(),gref->pyMode(),
2009  gref->pzMode(),gref->pMode());
2010  cand.setP4(p4);
2011  cand.setPositionAtECALEntrance(gsf->positionAtECALEntrance());
2012  } else if ( cfg_.produceEGCandsWithNoSuperCluster &&
2013  RO.primaryKFs.size() ) {
2014  const PFKFElement* kf = RO.primaryKFs[0].first;
2015  reco::TrackRef kref = RO.primaryKFs[0].first->trackRef();
2016  math::XYZTLorentzVector p4(kref->px(),kref->py(),kref->pz(),kref->p());
2017  cand.setP4(p4);
2018  cand.setPositionAtECALEntrance(kf->positionAtECALEntrance());
2019  }
2020  const float ele_mva_value = calculate_ele_mva(RO,xtra);
2021  fill_extra_info(RO,xtra);
2022  //std::cout << "PFEG ele_mva: " << ele_mva_value << std::endl;
2023  xtra.setMVA(ele_mva_value);
2024  cand.set_mva_e_pi(ele_mva_value);
2025  egcands.push_back(cand);
2026  egxs.push_back(xtra);
2027  }
2028 }
2029 
2030 float PFEGammaAlgo::
2031 calculate_ele_mva(const PFEGammaAlgo::ProtoEGObject& RO,
2032  reco::PFCandidateEGammaExtra& xtra) {
2033  if( !RO.primaryGSFs.size() ) return -2.0f;
2034  const PFGSFElement* gsfElement = RO.primaryGSFs.front().first;
2035  const PFKFElement* kfElement = NULL;
2036  if( RO.primaryKFs.size() ) kfElement = RO.primaryKFs.front().first;
2037  reco::GsfTrackRef RefGSF= gsfElement->GsftrackRef();
2038  reco::TrackRef RefKF;
2039  constexpr float m_el = 0.000511;
2040  const double Ein_gsf = std::hypot(RefGSF->pMode(),m_el);
2041  double deta_gsfecal = 1e6;
2042  double sigmaEtaEta = 1e-14;
2043  const double Ene_hcalgsf = std::accumulate(RO.hcalClusters.begin(),
2044  RO.hcalClusters.end(),
2045  0.0,
2046  [](const double a,
2047  const PFClusterFlaggedElement& b)
2048  { return a + b.first->clusterRef()->energy(); }
2049  );
2050  if( RO.primaryKFs.size() ) {
2051  RefKF = RO.primaryKFs.front().first->trackRef();
2052  }
2053  const double Eout_gsf = gsfElement->Pout().t();
2054  const double Etaout_gsf = gsfElement->positionAtECALEntrance().eta();
2055  double FirstEcalGsfEnergy(0.0), OtherEcalGsfEnergy(0.0), EcalBremEnergy(0.0);
2056  //shower shape of cluster closest to gsf track
2057  std::vector<const reco::PFCluster*> gsfcluster;
2058  for( const auto& ecal : RO.ecalclusters ) {
2059  const double cenergy = ecal.first->clusterRef()->correctedEnergy();
2060  ElementMap::value_type gsfToEcal(gsfElement,ecal.first);
2061  ElementMap::value_type kfToEcal(kfElement,ecal.first);
2062  bool hasgsf =
2063  ( std::find(RO.localMap.begin(), RO.localMap.end(), gsfToEcal) ==
2064  RO.localMap.end() );
2065  bool haskf =
2066  ( std::find(RO.localMap.begin(), RO.localMap.end(), kfToEcal) ==
2067  RO.localMap.end() );
2068  bool hasbrem = false;
2069  for( const auto& brem : RO.brems ) {
2070  ElementMap::value_type bremToEcal(brem.first,ecal.first);
2071  if( std::find(RO.localMap.begin(), RO.localMap.end(), bremToEcal) !=
2072  RO.localMap.end() ) {
2073  hasbrem = true;
2074  }
2075  }
2076  if( hasbrem && ecal.first != RO.electronClusters[0] ) {
2077  EcalBremEnergy += cenergy;
2078  }
2079  if( !hasbrem && ecal.first != RO.electronClusters[0] ) {
2080  if( hasgsf ) OtherEcalGsfEnergy += cenergy;
2081  if( haskf ) EcalBremEnergy += cenergy; // from conv. brem!
2082  if( !(hasgsf || haskf) ) OtherEcalGsfEnergy += cenergy; // stuff from SC
2083  }
2084  }
2085 
2086  if( RO.electronClusters[0] ) {
2087  reco::PFClusterRef cref = RO.electronClusters[0]->clusterRef();
2088  xtra.setGsfElectronClusterRef(_currentblock,*(RO.electronClusters[0]));
2089  FirstEcalGsfEnergy = cref->correctedEnergy();
2090  deta_gsfecal = cref->positionREP().eta() - Etaout_gsf;
2091  gsfcluster.push_back(&*cref);
2092  PFClusterWidthAlgo pfwidth(gsfcluster);
2093  sigmaEtaEta = pfwidth.pflowSigmaEtaEta();
2094  }
2095 
2096  // brem sequence information
2097  lateBrem = firstBrem = earlyBrem = -1.0f;
2098  if(RO.nBremsWithClusters > 0) {
2099  if (RO.lateBrem == 1) lateBrem = 1.0f;
2100  else lateBrem = 0.0f;
2101  firstBrem = RO.firstBrem;
2102  if(RO.firstBrem < 4) earlyBrem = 1.0f;
2103  else earlyBrem = 0.0f;
2104  }
2105  xtra.setEarlyBrem(earlyBrem);
2106  xtra.setLateBrem(lateBrem);
2107  if( FirstEcalGsfEnergy > 0.0 ) {
2108  if( RefGSF.isNonnull() ) {
2109  xtra.setGsfTrackPout(gsfElement->Pout());
2110  // normalization observables
2111  const float Pt_gsf = RefGSF->ptMode();
2112  lnPt_gsf = std::log(Pt_gsf);
2113  Eta_gsf = RefGSF->etaMode();
2114  // tracking observables
2115  const double ptModeErrorGsf = RefGSF->ptModeError();
2116  dPtOverPt_gsf = (ptModeErrorGsf > 0. ? ptModeErrorGsf/Pt_gsf : 1.0);
2117  nhit_gsf = RefGSF->hitPattern().trackerLayersWithMeasurement();
2118  chi2_gsf = RefGSF->normalizedChi2();
2119  DPtOverPt_gsf = (Pt_gsf - gsfElement->Pout().pt())/Pt_gsf;
2120  // kalman filter vars
2121  nhit_kf = 0;
2122  chi2_kf = -0.01;
2123  DPtOverPt_kf = -0.01;
2124  if( RefKF.isNonnull() ) {
2125  nhit_kf = RefKF->hitPattern().trackerLayersWithMeasurement();
2126  chi2_kf = RefKF->normalizedChi2();
2127  // not used for moment, weird behavior of variable
2128  // DPtOverPt_kf = (RefKF->pt() - RefKF->outerPt())/RefKF->pt();
2129  }
2130  //tracker + calorimetry observables
2131  const double EcalETot =
2132  (FirstEcalGsfEnergy+OtherEcalGsfEnergy+EcalBremEnergy);
2133  EtotPinMode = EcalETot / Ein_gsf;
2134  EGsfPoutMode = FirstEcalGsfEnergy / Eout_gsf;
2135  EtotBremPinPoutMode = ( (EcalBremEnergy + OtherEcalGsfEnergy) /
2136  (Ein_gsf - Eout_gsf) );
2137  DEtaGsfEcalClust = std::abs(deta_gsfecal);
2138  SigmaEtaEta = std::log(sigmaEtaEta);
2139  xtra.setDeltaEta(DEtaGsfEcalClust);
2140  xtra.setSigmaEtaEta(sigmaEtaEta);
2141 
2142  HOverHE = Ene_hcalgsf/(Ene_hcalgsf + FirstEcalGsfEnergy);
2143  HOverPin = Ene_hcalgsf / Ein_gsf;
2144  xtra.setHadEnergy(Ene_hcalgsf);
2145 
2146  // Apply bounds to variables and calculate MVA
2147  DPtOverPt_gsf = std::max(DPtOverPt_gsf,-0.2f);
2148  DPtOverPt_gsf = std::min(DPtOverPt_gsf,1.0f);
2149  dPtOverPt_gsf = std::min(dPtOverPt_gsf,0.3f);
2150  chi2_gsf = std::min(chi2_gsf,10.0f);
2151  DPtOverPt_kf = std::max(DPtOverPt_kf,-0.2f);
2152  DPtOverPt_kf = std::min(DPtOverPt_kf,1.0f);
2153  chi2_kf = std::min(chi2_kf,10.0f);
2154  EtotPinMode = std::max(EtotPinMode,0.0f);
2155  EtotPinMode = std::min(EtotPinMode,5.0f);
2156  EGsfPoutMode = std::max(EGsfPoutMode,0.0f);
2157  EGsfPoutMode = std::min(EGsfPoutMode,5.0f);
2158  EtotBremPinPoutMode = std::max(EtotBremPinPoutMode,0.0f);
2159  EtotBremPinPoutMode = std::min(EtotBremPinPoutMode,5.0f);
2160  DEtaGsfEcalClust = std::min(DEtaGsfEcalClust,0.1f);
2161  SigmaEtaEta = std::max(SigmaEtaEta,-14.0f);
2162  HOverPin = std::max(HOverPin,0.0f);
2163  HOverPin = std::min(HOverPin,5.0f);
2164  /*
2165  std::cout << " **** PFEG BDT observables ****" << endl;
2166  std::cout << " < Normalization > " << endl;
2167  std::cout << " Pt_gsf " << Pt_gsf << " Pin " << Ein_gsf
2168  << " Pout " << Eout_gsf << " Eta_gsf " << Eta_gsf << endl;
2169  std::cout << " < PureTracking > " << endl;
2170  std::cout << " dPtOverPt_gsf " << dPtOverPt_gsf
2171  << " DPtOverPt_gsf " << DPtOverPt_gsf
2172  << " chi2_gsf " << chi2_gsf
2173  << " nhit_gsf " << nhit_gsf
2174  << " DPtOverPt_kf " << DPtOverPt_kf
2175  << " chi2_kf " << chi2_kf
2176  << " nhit_kf " << nhit_kf << endl;
2177  std::cout << " < track-ecal-hcal-ps " << endl;
2178  std::cout << " EtotPinMode " << EtotPinMode
2179  << " EGsfPoutMode " << EGsfPoutMode
2180  << " EtotBremPinPoutMode " << EtotBremPinPoutMode
2181  << " DEtaGsfEcalClust " << DEtaGsfEcalClust
2182  << " SigmaEtaEta " << SigmaEtaEta
2183  << " HOverHE " << HOverHE << " Hcal energy " << Ene_hcalgsf
2184  << " HOverPin " << HOverPin
2185  << " lateBrem " << lateBrem
2186  << " firstBrem " << firstBrem << endl;
2187  */
2188 
2189  return tmvaReaderEle_->EvaluateMVA("BDT");
2190  }
2191  }
2192  return -2.0f;
2193 }
2194 
2195 void PFEGammaAlgo::fill_extra_info( const ProtoEGObject& RO,
2196  reco::PFCandidateEGammaExtra& xtra ) {
2197  // add tracks associated to clusters that are not T_FROM_GAMMACONV
2198  // info about single-leg convs is already save, so just veto in loops
2199  IsConversionTrack<reco::PFBlockElementTrack> isConvKf;
2200  auto KFbegin = _splayedblock[reco::PFBlockElement::TRACK].begin();
2201  auto KFend = _splayedblock[reco::PFBlockElement::TRACK].end();
2202  for( auto& ecal : RO.ecalclusters ) {
2203  NotCloserToOther<reco::PFBlockElement::ECAL,
2204  reco::PFBlockElement::TRACK,
2205  true>
2206  ECALToTracks(_currentblock,_currentlinks,ecal.first);
2207  auto notmatchedkf = std::partition(KFbegin,KFend,ECALToTracks);
2208  auto notconvkf = std::partition(KFbegin,notmatchedkf,isConvKf);
2209  // go through non-conv-identified kfs and check MVA to add conversions
2210  for( auto kf = notconvkf; kf != notmatchedkf; ++kf ) {
2211  const reco::PFBlockElementTrack* elemaskf =
2212  docast(const reco::PFBlockElementTrack*,kf->first);
2213  xtra.addExtraNonConvTrack(_currentblock,*elemaskf);
2214  }
2215  }
2216 }
2217 
2218 // currently stolen from PFECALSuperClusterAlgo, we should
2219 // try to factor this correctly since the operation is the same in
2220 // both places...
2221 reco::SuperCluster PFEGammaAlgo::
2222 buildRefinedSuperCluster(const PFEGammaAlgo::ProtoEGObject& RO) {
2223  if( !RO.ecalclusters.size() ) {
2224  return reco::SuperCluster(0.0,math::XYZPoint(0,0,0));
2225  }
2226 
2227  SumPSEnergy sumps1(reco::PFBlockElement::PS1),
2228  sumps2(reco::PFBlockElement::PS2);
2229 
2230  bool isEE = false;
2231  edm::Ptr<reco::PFCluster> clusptr;
2232  // need the vector of raw pointers for a PF width class
2233  std::vector<const reco::PFCluster*> bare_ptrs;
2234  // calculate necessary parameters and build the SC
2235  double posX(0), posY(0), posZ(0),
2236  rawSCEnergy(0), corrSCEnergy(0), corrPSEnergy(0),
2237  PS1_clus_sum(0), PS2_clus_sum(0),
2238  ePS1(0), ePS2(0), ps1_energy(0.0), ps2_energy(0.0);
2239  for( auto& clus : RO.ecalclusters ) {
2240  ePS1 = 0;
2241  ePS2 = 0;
2242  isEE = PFLayer::ECAL_ENDCAP == clus.first->clusterRef()->layer();
2243  clusptr =
2244  edm::refToPtr<reco::PFClusterCollection>(clus.first->clusterRef());
2245  bare_ptrs.push_back(clusptr.get());
2246 
2247  const double cluseraw = clusptr->energy();
2248  double cluscalibe = clusptr->correctedEnergy();
2249  const math::XYZPoint& cluspos = clusptr->position();
2250  posX += cluseraw * cluspos.X();
2251  posY += cluseraw * cluspos.Y();
2252  posZ += cluseraw * cluspos.Z();
2253  // update EE calibrated super cluster energies
2254  if( isEE ) {
2255  const auto& psclusters = RO.ecal2ps.at(clus.first);
2256  PS1_clus_sum = std::accumulate(psclusters.begin(),psclusters.end(),
2257  0.0,sumps1);
2258  PS2_clus_sum = std::accumulate(psclusters.begin(),psclusters.end(),
2259  0.0,sumps2);
2260  cluscalibe =
2261  cfg_.thePFEnergyCalibration->energyEm(*clusptr,
2262  PS1_clus_sum,PS2_clus_sum,
2263  ePS1, ePS2,
2264  cfg_.applyCrackCorrections);
2265  }
2266 
2267  rawSCEnergy += cluseraw;
2268  corrSCEnergy += cluscalibe;
2269  ps1_energy += ePS1;
2270  ps2_energy += ePS2;
2271  corrPSEnergy += ePS1 + ePS2;
2272  }
2273  posX /= rawSCEnergy;
2274  posY /= rawSCEnergy;
2275  posZ /= rawSCEnergy;
2276 
2277  // now build the supercluster
2278  reco::SuperCluster new_sc(corrSCEnergy,math::XYZPoint(posX,posY,posZ));
2279 
2280  clusptr =
2281  edm::refToPtr<reco::PFClusterCollection>(RO.ecalclusters.front().
2282  first->clusterRef());
2283  new_sc.setCorrectedEnergy(corrSCEnergy);
2284  new_sc.setSeed(clusptr);
2285  new_sc.setPreshowerEnergyPlane1(ps1_energy);
2286  new_sc.setPreshowerEnergyPlane2(ps2_energy);
2287  new_sc.setPreshowerEnergy(corrPSEnergy);
2288  for( const auto& clus : RO.ecalclusters ) {
2289  clusptr =
2290  edm::refToPtr<reco::PFClusterCollection>(clus.first->clusterRef());
2291  new_sc.addCluster(clusptr);
2292  auto& hits_and_fractions = clusptr->hitsAndFractions();
2293  for( auto& hit_and_fraction : hits_and_fractions ) {
2294  new_sc.addHitAndFraction(hit_and_fraction.first,hit_and_fraction.second);
2295  }
2296  // put the preshower stuff back in later
2297  const auto& cluspsassociation = RO.ecal2ps.at(clus.first);
2298  // EE rechits should be uniquely matched to sets of pre-shower
2299  // clusters at this point, so we throw an exception if otherwise
2300  // now wrapped in EDM debug flags
2301  for( const auto& pscluselem : cluspsassociation ) {
2302  edm::Ptr<reco::PFCluster> psclus =
2303  edm::refToPtr<reco::PFClusterCollection>(pscluselem.first->
2304  clusterRef());
2305 #ifdef PFFLOW_DEBUG
2306  auto found_pscluster = std::find(new_sc.preshowerClustersBegin(),
2307  new_sc.preshowerClustersEnd(),
2308  reco::CaloClusterPtr(psclus));
2309  if( found_pscluster == new_sc.preshowerClustersEnd() ) {
2310 #endif
2311  new_sc.addPreshowerCluster(psclus);
2312 #ifdef PFFLOW_DEBUG
2313  } else {
2314  throw cms::Exception("PFECALSuperClusterAlgo::buildSuperCluster")
2315  << "Found a PS cluster matched to more than one EE cluster!"
2316  << std::endl << std::hex << psclus.get() << " == "
2317  << found_pscluster->get() << std::dec << std::endl;
2318  }
2319 #endif
2320  }
2321  }
2322 
2323  // calculate linearly weighted cluster widths
2324  PFClusterWidthAlgo pfwidth(bare_ptrs);
2325  new_sc.setEtaWidth(pfwidth.pflowEtaWidth());
2326  new_sc.setPhiWidth(pfwidth.pflowPhiWidth());
2327 
2328  // cache the value of the raw energy
2329  new_sc.rawEnergy();
2330 
2331  return new_sc;
2332 }
2333 
2334 void PFEGammaAlgo::
2335 unlinkRefinableObjectKFandECALWithBadEoverP(ProtoEGObject& RO) {
2336  // this only means something for ROs with a primary GSF track
2337  if( !RO.primaryGSFs.size() ) return;
2338  // need energy sums to tell if we've added crap or not
2339  const double Pin_gsf = RO.primaryGSFs.front().first->GsftrackRef()->pMode();
2340  const double gsfOuterEta =
2341  RO.primaryGSFs.front().first->positionAtECALEntrance().Eta();
2342  double tot_ecal= 0.0;
2343  std::vector<double> min_brem_dists;
2344  std::vector<double> closest_brem_eta;
2345  // first get the total ecal energy (we should replace this with a cache)
2346  for( const auto& ecal : RO.ecalclusters ) {
2347  tot_ecal += ecal.first->clusterRef()->correctedEnergy();
2348  // we also need to look at the minimum distance to brems
2349  // since energetic brems will be closer to the brem than the track
2350  double min_brem_dist = 5000.0;
2351  double eta = -999.0;
2352  for( const auto& brem : RO.brems ) {
2353  const float dist = _currentblock->dist(brem.first->index(),
2354  ecal.first->index(),
2355  _currentlinks,
2356  reco::PFBlock::LINKTEST_ALL);
2357  if( dist < min_brem_dist && dist != -1.0f ) {
2358  min_brem_dist = dist;
2359  eta = brem.first->positionAtECALEntrance().Eta();
2360  }
2361  }
2362  min_brem_dists.push_back(min_brem_dist);
2363  closest_brem_eta.push_back(eta);
2364  }
2365 
2366  // loop through the ECAL clusters and remove ECAL clusters matched to
2367  // secondary track either in *or* out of the SC if the E/pin is bad
2368  for( auto secd_kf = RO.secondaryKFs.begin();
2369  secd_kf != RO.secondaryKFs.end(); ++secd_kf ) {
2370  reco::TrackRef trkRef = secd_kf->first->trackRef();
2371  const float secpin = secd_kf->first->trackRef()->p();
2372  bool remove_this_kf = false;
2373  for( auto ecal = RO.ecalclusters.begin();
2374  ecal != RO.ecalclusters.end(); ++ecal ) {
2375  size_t bremidx = std::distance(RO.ecalclusters.begin(),ecal);
2376  const float minbremdist = min_brem_dists[bremidx];
2377  const double ecalenergy = ecal->first->clusterRef()->correctedEnergy();
2378  const double Epin = ecalenergy/secpin;
2379  const double detaGsf =
2380  std::abs(gsfOuterEta - ecal->first->clusterRef()->positionREP().Eta());
2381  const double detaBrem =
2382  std::abs(closest_brem_eta[bremidx] -
2383  ecal->first->clusterRef()->positionREP().Eta());
2384 
2385  ElementMap::value_type check_match(ecal->first,secd_kf->first);
2386  auto kf_matched = std::find(RO.localMap.begin(),
2387  RO.localMap.end(),
2388  check_match);
2389 
2390  const float tkdist = _currentblock->dist(secd_kf->first->index(),
2391  ecal->first->index(),
2392  _currentlinks,
2393  reco::PFBlock::LINKTEST_ALL);
2394 
2395  // do not reject this track if it is closer to a brem than the
2396  // secondary track, or if it lies in the delta-eta plane with the
2397  // gsf track or if it is in the dEta plane with the brems
2398  if( Epin > 3 && kf_matched != RO.localMap.end() &&
2399  tkdist != -1.0f && tkdist < minbremdist &&
2400  detaGsf > 0.05 && detaBrem > 0.015) {
2401  double res_with = std::abs((tot_ecal-Pin_gsf)/Pin_gsf);
2402  double res_without = std::abs((tot_ecal-ecalenergy-Pin_gsf)/Pin_gsf);
2403  if(res_without < res_with) {
2404  LOGDRESSED("PFEGammaAlgo")
2405  << " REJECTED_RES totenergy " << tot_ecal
2406  << " Pin_gsf " << Pin_gsf
2407  << " cluster to secondary " << ecalenergy
2408  << " res_with " << res_with
2409  << " res_without " << res_without << std::endl;
2410  tot_ecal -= ecalenergy;
2411  remove_this_kf = true;
2412  ecal = RO.ecalclusters.erase(ecal);
2413  if( ecal == RO.ecalclusters.end() ) break;
2414  }
2415  }
2416  }
2417  if( remove_this_kf ) {
2418  secd_kf = RO.secondaryKFs.erase(secd_kf);
2419  if( secd_kf == RO.secondaryKFs.end() ) break;
2420  }
2421  }
2422 }
2423 
2424 void PFEGammaAlgo::
2425 unlinkRefinableObjectKFandECALMatchedToHCAL(ProtoEGObject& RO,
2426  bool removeFreeECAL,
2427  bool removeSCEcal) {
2428  std::vector<bool> cluster_in_sc;
2429  auto ecal_begin = RO.ecalclusters.begin();
2430  auto ecal_end = RO.ecalclusters.end();
2431  auto hcal_begin = _splayedblock[reco::PFBlockElement::HCAL].begin();
2432  auto hcal_end = _splayedblock[reco::PFBlockElement::HCAL].end();
2433  for( auto secd_kf = RO.secondaryKFs.begin();
2434  secd_kf != RO.secondaryKFs.end(); ++secd_kf ) {
2435  bool remove_this_kf = false;
2436  NotCloserToOther<reco::PFBlockElement::TRACK,reco::PFBlockElement::HCAL>
2437  tracksToHCALs(_currentblock,_currentlinks,secd_kf->first);
2438  reco::TrackRef trkRef = secd_kf->first->trackRef();
2439  const unsigned int Algo = whichTrackAlgo(trkRef);
2440  const float secpin = trkRef->p();
2441 
2442  for( auto ecal = ecal_begin; ecal != ecal_end; ++ecal ) {
2443  const double ecalenergy = ecal->first->clusterRef()->correctedEnergy();
2444  // first check if the cluster is in the SC (use dist calc for fastness)
2445  const size_t clus_idx = std::distance(ecal_begin,ecal);
2446  if( cluster_in_sc.size() < clus_idx + 1) {
2447  float dist = -1.0f;
2448  if( RO.parentSC ) {
2449  dist = _currentblock->dist(secd_kf->first->index(),
2450  ecal->first->index(),
2451  _currentlinks,
2452  reco::PFBlock::LINKTEST_ALL);
2453  }
2454  cluster_in_sc.push_back(dist != -1.0f);
2455  }
2456 
2457  ElementMap::value_type check_match(ecal->first,secd_kf->first);
2458  auto kf_matched = std::find(RO.localMap.begin(),
2459  RO.localMap.end(),
2460  check_match);
2461  // if we've found a secondary KF that matches this ecal cluster
2462  // now we see if it is matched to HCAL
2463  // if it is matched to an HCAL cluster we take different
2464  // actions if the cluster was in an SC or not
2465  if( kf_matched != RO.localMap.end() ) {
2466  auto hcal_matched = std::partition(hcal_begin,hcal_end,tracksToHCALs);
2467  for( auto hcalclus = hcal_begin;
2468  hcalclus != hcal_matched;
2469  ++hcalclus ) {
2470  const reco::PFBlockElementCluster * clusthcal =
2471  dynamic_cast<const reco::PFBlockElementCluster*>(hcalclus->first);
2472  const double hcalenergy = clusthcal->clusterRef()->energy();
2473  const double hpluse = ecalenergy+hcalenergy;
2474  const bool isHoHE = ( (hcalenergy / hpluse ) > 0.1 && Algo < 3 );
2475  const bool isHoE = ( hcalenergy > ecalenergy );
2476  const bool isPoHE = ( secpin > hpluse );
2477  if( cluster_in_sc[clus_idx] ) {
2478  if(isHoE || isPoHE) {
2479  LOGDRESSED("PFEGammaAlgo")
2480  << "REJECTED TRACK FOR H/E or P/(H+E), CLUSTER IN SC"
2481  << " H/H+E " << (hcalenergy / hpluse)
2482  << " H/E " << (hcalenergy > ecalenergy)
2483  << " P/(H+E) " << (secpin/hpluse)
2484  << " HCAL ENE " << hcalenergy
2485  << " ECAL ENE " << ecalenergy
2486  << " secPIN " << secpin
2487  << " Algo Track " << Algo << std::endl;
2488  remove_this_kf = true;
2489  }
2490  } else {
2491  if(isHoHE){
2492  LOGDRESSED("PFEGammaAlgo")
2493  << "REJECTED TRACK FOR H/H+E, CLUSTER NOT IN SC"
2494  << " H/H+E " << (hcalenergy / hpluse)
2495  << " H/E " << (hcalenergy > ecalenergy)
2496  << " P/(H+E) " << (secpin/hpluse)
2497  << " HCAL ENE " << hcalenergy
2498  << " ECAL ENE " << ecalenergy
2499  << " secPIN " << secpin
2500  << " Algo Track " << Algo << std::endl;
2501  remove_this_kf = true;
2502  }
2503  }
2504  }
2505  }
2506  }
2507  if( remove_this_kf ) {
2508  secd_kf = RO.secondaryKFs.erase(secd_kf);
2509  if( secd_kf == RO.secondaryKFs.end() ) break;
2510  }
2511  }
2512 }
2513 
2514 
2515 unsigned int PFEGammaAlgo::whichTrackAlgo(const reco::TrackRef& trackRef) {
2516  unsigned int Algo = 0;
2517  switch (trackRef->algo()) {
2518  case TrackBase::ctf:
2519  case TrackBase::iter0:
2520  case TrackBase::iter1:
2521  case TrackBase::iter2:
2522  Algo = 0;
2523  break;
2524  case TrackBase::iter3:
2525  Algo = 1;
2526  break;
2527  case TrackBase::iter4:
2528  Algo = 2;
2529  break;
2530  case TrackBase::iter5:
2531  Algo = 3;
2532  break;
2533  case TrackBase::iter6:
2534  Algo = 4;
2535  break;
2536  default:
2537  Algo = 5;
2538  break;
2539  }
2540  return Algo;
2541 }
2542 bool PFEGammaAlgo::isPrimaryTrack(const reco::PFBlockElementTrack& KfEl,
2543  const reco::PFBlockElementGsfTrack& GsfEl) {
2544  bool isPrimary = false;
2545 
2546  GsfPFRecTrackRef gsfPfRef = GsfEl.GsftrackRefPF();
2547 
2548  if(gsfPfRef.isNonnull()) {
2549  PFRecTrackRef kfPfRef = KfEl.trackRefPF();
2550  PFRecTrackRef kfPfRef_fromGsf = (*gsfPfRef).kfPFRecTrackRef();
2551  if(kfPfRef.isNonnull() && kfPfRef_fromGsf.isNonnull()) {
2552  reco::TrackRef kfref= (*kfPfRef).trackRef();
2553  reco::TrackRef kfref_fromGsf = (*kfPfRef_fromGsf).trackRef();
2554  if(kfref.isNonnull() && kfref_fromGsf.isNonnull()) {
2555  if(kfref == kfref_fromGsf)
2556  isPrimary = true;
2557  }
2558  }
2559  }
2560 
2561  return isPrimary;
2562 }
PFEGammaAlgo::EtotBremPinPoutMode
float EtotBremPinPoutMode
Definition: PFEGammaAlgo.h:309
type
type
Definition: HCALResponse.h:21
reco::PFBlockElementSuperCluster::superClusterRef
const SuperClusterRef & superClusterRef() const
Definition: PFBlockElementSuperCluster.h:37
reco::PFBlockElementGsfTrack::Pout
const math::XYZTLorentzVector & Pout() const
Definition: PFBlockElementGsfTrack.h:70
PFEGammaAlgo::mvaValue
double mvaValue
Definition: PFEGammaAlgo.h:350
reco::PFBlockElement
Abstract base class for a PFBlock element (track, cluster...)
Definition: PFBlockElement.h:28
reco::CaloCluster::position
const math::XYZPoint & position() const
cluster centroid position
Definition: CaloCluster.h:124
reco::SuperCluster::preshowerClustersBegin
CaloCluster_iterator preshowerClustersBegin() const
fist iterator over PreshowerCluster constituents
Definition: SuperCluster.h:81
i
int i
Definition: DBlmapReader.cc:9
reco::PFBlockElementTrack::trackRef
const reco::TrackRef & trackRef() const
Definition: PFBlockElementTrack.h:57
PFEGammaAlgo::whichTrackAlgo
unsigned int whichTrackAlgo(const reco::TrackRef &trackRef)
Definition: PFEGammaAlgo.cc:2515
reco::LeafCandidate::setCharge
virtual void setCharge(Charge q) GCC11_FINAL
set electric charge
Definition: LeafCandidate.h:142
reco::PFBlockElementGsfTrack::GsftrackRef
const reco::GsfTrackRef & GsftrackRef() const
Definition: PFBlockElementGsfTrack.h:53
info
static const TGPicture * info(bool iBackgroundIsBlack)
Definition: FWCollectionSummaryWidget.cc:170
PFEGammaAlgo::ProtoEGObject::ecalclusters
std::vector< PFClusterFlaggedElement > ecalclusters
Definition: PFEGammaAlgo.h:83
ClusterClusterMapping::overlap
static bool overlap(const reco::CaloCluster &sc1, const reco::CaloCluster &sc, float minfrac=0.01, bool debug=false)
Definition: ClusterClusterMapping.cc:4
reco::PFBlockElementGsfTrack::Dump
void Dump(std::ostream &out=std::cout, const char *tab=" ") const
print the object inside the element
Definition: PFBlockElementGsfTrack.cc:37
reco::PFCandidateEGammaExtra::setSuperClusterRef
void setSuperClusterRef(reco::SuperClusterRef sc)
set reference to the corresponding supercluster
Definition: PFCandidateEGammaExtra.h:111
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
dbtoconf.parent
list parent
Definition: dbtoconf.py:74
PFEGammaAlgo::egExtra_
reco::PFCandidateEGammaExtraCollection egExtra_
Definition: PFEGammaAlgo.h:384
PFEGammaAlgo::fifthStepKfTrack_
std::vector< std::pair< unsigned int, unsigned int > > fifthStepKfTrack_
Definition: PFEGammaAlgo.h:289
PFEGammaAlgo::refinedscs_
reco::SuperClusterCollection refinedscs_
Definition: PFEGammaAlgo.h:178
reco::CaloCluster::addHitAndFraction
void addHitAndFraction(DetId id, float fraction)
Definition: CaloCluster.h:183
PFEGammaAlgo::RunPFEG
void RunPFEG(const reco::PFBlockRef &blockRef, std::vector< bool > &active)
Definition: PFEGammaAlgo.cc:755
PFEGammaAlgo::ProtoEGObject::parentBlock
reco::PFBlockRef parentBlock
Definition: PFEGammaAlgo.h:77
PFEGammaAlgo::unlinkRefinableObjectKFandECALWithBadEoverP
void unlinkRefinableObjectKFandECALWithBadEoverP(ProtoEGObject &)
Definition: PFEGammaAlgo.cc:2335
reco::Vertex::tracks_end
trackRef_iterator tracks_end() const
last iterator over tracks
Definition: Vertex.cc:44
PFMuonAlgo::isMuon
static bool isMuon(const reco::PFBlockElement &elt)
Definition: PFMuonAlgo.cc:153
reco::PFCandidate::setPositionAtECALEntrance
void setPositionAtECALEntrance(const math::XYZPointF &pos)
set position at ECAL entrance
Definition: PFCandidate.h:330
reco::PFCandidateEGammaExtra::setGsfElectronClusterRef
void setGsfElectronClusterRef(const reco::PFBlockRef &blk, const reco::PFBlockElementCluster &ref)
set gsf electron cluster ref
Definition: PFCandidateEGammaExtra.h:88
PFEGammaAlgo::EGsfPoutMode
float EGsfPoutMode
Definition: PFEGammaAlgo.h:309
reco::PFCandidateEGammaExtra::addSingleLegConvTrackRefMva
void addSingleLegConvTrackRefMva(const std::pair< reco::TrackRef, float > &trackrefmva)
add Single Leg Conversion TrackRef
Definition: PFCandidateEGammaExtra.cc:156
reco::SuperCluster::setPreshowerEnergyPlane2
void setPreshowerEnergyPlane2(double preshowerEnergy2)
Definition: SuperCluster.h:61
PFEGammaAlgo::PFClusterFlaggedElement
std::pair< const PFClusterElement *, bool > PFClusterFlaggedElement
Definition: PFEGammaAlgo.h:64
reco::Vertex::y
double y() const
y coordinate
Definition: Vertex.h:110
reco::PFCandidateEGammaExtra::setHadEnergy
void setHadEnergy(float val)
set the had energy. The cluster energies should be entered before
Definition: PFCandidateEGammaExtra.cc:79
reco::PFBlockElementGsfTrack::positionAtECALEntrance
const math::XYZPointF & positionAtECALEntrance() const
Definition: PFBlockElementGsfTrack.h:61
PFEGammaAlgo::PFEGConfigInfo::thePFEnergyCalibration
std::shared_ptr< PFEnergyCalibration > thePFEnergyCalibration
Definition: PFEGammaAlgo.h:107
reco::PFBlockElement::type
Type type() const
Definition: PFBlockElement.h:83
reco::PFBlock::LinkData
std::map< unsigned int, Link > LinkData
Definition: PFBlock.h:46
PV3DBase::phi
Geom::Phi< T > phi() const
Definition: PV3DBase.h:69
PFClusterWidthAlgo::pflowPhiWidth
double pflowPhiWidth() const
Definition: PFClusterWidthAlgo.h:16
PFEGammaAlgo::ProtoEGObject::hcalClusters
std::vector< PFClusterFlaggedElement > hcalClusters
Definition: PFEGammaAlgo.h:95
PFEGammaAlgo::isPrimaryTrack
bool isPrimaryTrack(const reco::PFBlockElementTrack &KfEl, const reco::PFBlockElementGsfTrack &GsfEl)
Definition: PFEGammaAlgo.cc:2542
PFEGammaAlgo::DEtaGsfEcalClust
float DEtaGsfEcalClust
Definition: PFEGammaAlgo.h:310
X
#define X(str)
Definition: MuonsGrabber.cc:48
fff_deleter.log
tuple log
Definition: fff_deleter.py:13
PFEGammaAlgo::linkRefinableObjectPrimaryGSFTrackToECAL
void linkRefinableObjectPrimaryGSFTrackToECAL(ProtoEGObject &)
Definition: PFEGammaAlgo.cc:1596
asciidump.elements
list elements
Definition: asciidump.py:414
NULL
#define NULL
Definition: scimark2.h:8
reco::PFBlockElementGsfTrack::GsftrackRefPF
const GsfPFRecTrackRef & GsftrackRefPF() const
Definition: PFBlockElementGsfTrack.h:48
reco::PFBlockElementCluster::clusterRef
const PFClusterRef & clusterRef() const
Definition: PFBlockElementCluster.h:32
reco::PFBlock::elements
const edm::OwnVector< reco::PFBlockElement > & elements() const
Definition: PFBlock.h:107
Phi_mpi_pi
double Phi_mpi_pi(double x)
Definition: JetUtil.h:24
reco::TrackExtraRef
edm::Ref< TrackExtraCollection > TrackExtraRef
persistent reference to a TrackExtra
Definition: TrackExtraFwd.h:13
PFEGammaAlgo::ProtoEGObject::electronClusters
std::vector< const PFClusterElement * > electronClusters
Definition: PFEGammaAlgo.h:99
reco::Vertex::index
unsigned int index
index type
Definition: Vertex.h:49
reco::PFBlock::linkData
const LinkData & linkData() const
Definition: PFBlock.h:112
PFEGammaAlgo::_splayedblock
std::vector< std::vector< PFFlaggedElement > > _splayedblock
Definition: PFEGammaAlgo.h:187
reco::SuperCluster::setSeed
void setSeed(const CaloClusterPtr &r)
list of used xtals by DetId // now inherited by CaloCluster
Definition: SuperCluster.h:96
reco::PFBlockElementGsfTrack::Pin
const math::XYZTLorentzVector & Pin() const
Definition: PFBlockElementGsfTrack.h:69
PFEGammaAlgo::PFGSFFlaggedElement
std::pair< const PFGSFElement *, bool > PFGSFFlaggedElement
Definition: PFEGammaAlgo.h:62
reco::CaloClusterPtr
edm::Ptr< CaloCluster > CaloClusterPtr
Definition: CaloClusterFwd.h:17
PFEGammaAlgo::linkRefinableObjectBremTangentsToECAL
void linkRefinableObjectBremTangentsToECAL(ProtoEGObject &)
Definition: PFEGammaAlgo.cc:1739
math::XYZPointF
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< float > > XYZPointF
point in space with cartesian internal representation
Definition: Point3D.h:10
edm::Ref::isAvailable
bool isAvailable() const
Definition: Ref.h:276
PFEGammaAlgo::X0_outer
TH2D * X0_outer
Definition: PFEGammaAlgo.h:373
PFEGammaAlgo::PFEGConfigInfo::primaryVtx
const reco::Vertex * primaryVtx
Definition: PFEGammaAlgo.h:124
eta
T eta() const
Definition: Basic3DVectorLD.h:177
reco::PFBlockElementBrem::indTrajPoint
unsigned int indTrajPoint() const
Definition: PFBlockElementBrem.h:45
AlCaHLTBitMon_QueryRunRegistry.comp
string comp
Definition: AlCaHLTBitMon_QueryRunRegistry.py:248
dt_dqm_sourceclient_common_cff.reco
tuple reco
Definition: dt_dqm_sourceclient_common_cff.py:107
spr::find
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:7
PFEGammaAlgo::_currentblock
reco::PFBlockRef _currentblock
Definition: PFEGammaAlgo.h:183
reco::PFCandidateEGammaExtra::setGsfTrackRef
void setGsfTrackRef(const reco::GsfTrackRef &ref)
set gsftrack reference
Definition: PFCandidateEGammaExtra.cc:47
PFEGammaAlgo::PFFlaggedElement
std::pair< const reco::PFBlockElement *, bool > PFFlaggedElement
Definition: PFEGammaAlgo.h:59
reco::SuperCluster::setPhiWidth
void setPhiWidth(double pw)
Definition: SuperCluster.h:62
reco::Vertex::position
const Point & position() const
position
Definition: Vertex.h:106
PFClusterWidthAlgo::pflowEtaWidth
double pflowEtaWidth() const
Definition: PFClusterWidthAlgo.h:17
PFEGammaAlgo::PFSCFlaggedElement
std::pair< const PFSCElement *, bool > PFSCFlaggedElement
Definition: PFEGammaAlgo.h:60
reco::PFBlockElementTrack::positionAtECALEntrance
const math::XYZPointF & positionAtECALEntrance() const
Definition: PFBlockElementTrack.h:48
PFEGammaAlgo::cfg_
PFEGConfigInfo cfg_
Definition: PFEGammaAlgo.h:292
reco::PFCandidate::set_mva_e_pi
void set_mva_e_pi(float mva)
Definition: PFCandidate.h:291
PFEGammaAlgo::DPtOverPt_gsf
float DPtOverPt_gsf
Definition: PFEGammaAlgo.h:303
PFEGammaAlgo::linkRefinableObjectECALToSingleLegConv
void linkRefinableObjectECALToSingleLegConv(ProtoEGObject &)
Definition: PFEGammaAlgo.cc:1842
PFEGammaAlgo::eetops_
edm::Handle< reco::PFCluster::EEtoPSAssociation > eetops_
Definition: PFEGammaAlgo.h:182
math::XYZTLorentzVector
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:29
PFEGammaAlgo::EvaluateSingleLegMVA
float EvaluateSingleLegMVA(const reco::PFBlockRef &blockref, const reco::Vertex &primaryvtx, unsigned int track_index)
Definition: PFEGammaAlgo.cc:771
edm::Ref::isNonnull
bool isNonnull() const
Checks for non-null.
Definition: Ref.h:250
LOGERR
#define LOGERR(x)
Definition: PFEGammaAlgo.cc:45
PFEGammaAlgo::ProtoEGObject::secondaryKFs
std::vector< PFKFFlaggedElement > secondaryKFs
Definition: PFEGammaAlgo.h:92
reco::SuperCluster::setEtaWidth
void setEtaWidth(double ew)
Definition: SuperCluster.h:63
reco::PFCandidateEGammaExtra::setSigmaEtaEta
void setSigmaEtaEta(float val)
set the sigmaetaeta
Definition: PFCandidateEGammaExtra.cc:85
LOGVERB
#define LOGVERB(x)
Definition: PFEGammaAlgo.cc:43
LOGDRESSED
#define LOGDRESSED(x)
Definition: PFEGammaAlgo.cc:46
PFEGammaAlgo::dumpCurrentRefinableObjects
void dumpCurrentRefinableObjects() const
Definition: PFEGammaAlgo.cc:1285
edm::Ref::isNull
bool isNull() const
Checks for null.
Definition: Ref.h:247
PFEGammaAlgo::initializeProtoCands
void initializeProtoCands(std::list< ProtoEGObject > &)
Definition: PFEGammaAlgo.cc:979
edm::Ptr::get
T const * get() const
Returns C++ pointer to the item.
Definition: Ptr.h:143
PFEGammaAlgo::linkRefinableObjectPrimaryKFsToSecondaryKFs
void linkRefinableObjectPrimaryKFsToSecondaryKFs(ProtoEGObject &)
Definition: PFEGammaAlgo.cc:1560
reco::PFBlockElement::index
unsigned index() const
Definition: PFBlockElement.h:99
max
const T & max(const T &a, const T &b)
Definition: MaterialBudgetTrackerHistos.cc:4
PFEGammaAlgo::unwrapSuperCluster
bool unwrapSuperCluster(const reco::PFBlockElementSuperCluster *, std::vector< PFClusterFlaggedElement > &, ClusterMap &)
Definition: PFEGammaAlgo.cc:1132
PFClusterWidthAlgo::pflowSigmaEtaEta
double pflowSigmaEtaEta() const
Definition: PFClusterWidthAlgo.h:18
reco::CaloCluster::setCorrectedEnergy
void setCorrectedEnergy(double cenergy)
Definition: CaloCluster.h:109
reco::PFCandidate::addElementInBlock
void addElementInBlock(const reco::PFBlockRef &blockref, unsigned elementIndex)
add an element to the current PFCandidate
Definition: PFCandidate.cc:207
std::swap
void swap(edm::DataFrameContainer &lhs, edm::DataFrameContainer &rhs)
Definition: DataFrameContainer.h:223
PFEGammaAlgo::EEtoPSAssociation
reco::PFCluster::EEtoPSAssociation EEtoPSAssociation
Definition: PFEGammaAlgo.h:53
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:48
p4
double p4[4]
Definition: TauolaWrapper.h:92
reco::PFBlockElementTrack::convRefs
const ConversionRefVector & convRefs() const
Definition: PFBlockElementTrack.h:103
PFEGammaAlgo::PFEGammaAlgo
PFEGammaAlgo(const PFEGConfigInfo &)
Definition: PFEGammaAlgo.cc:687
PFEGammaAlgo::ProtoEGObject::parentSC
const PFSCElement * parentSC
Definition: PFEGammaAlgo.h:78
query.result
tuple result
Definition: query.py:137
reco::PFCandidateEGammaExtraCollection
std::vector< reco::PFCandidateEGammaExtra > PFCandidateEGammaExtraCollection
collection of PFCandidateEGammaExtras
Definition: PFCandidateEGammaExtraFwd.h:12
PFEGammaAlgo::ProtoEGObject::brems
std::vector< PFBremFlaggedElement > brems
Definition: PFEGammaAlgo.h:89
reco::PFCandidateEGammaExtra::setEarlyBrem
void setEarlyBrem(float val)
set EarlyBrem
Definition: PFCandidateEGammaExtra.cc:74
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
reco::Vertex::z
double z() const
y coordinate
Definition: Vertex.h:112
PFEGammaAlgo::DPtOverPt_kf
float DPtOverPt_kf
Definition: PFEGammaAlgo.h:303
reco::PFCandidateEGammaExtra::setGsfTrackPout
void setGsfTrackPout(const math::XYZTLorentzVector &pout)
set the pout (not trivial to get from the GSF track)
Definition: PFCandidateEGammaExtra.cc:51
PFEGammaAlgo::dPtOverPt_gsf
float dPtOverPt_gsf
Definition: PFEGammaAlgo.h:303
f
double f[11][100]
Definition: MuScleFitUtils.cc:78
reco::CaloCluster::energy
double energy() const
cluster energy
Definition: CaloCluster.h:120
end
#define end
Definition: vmac.h:37
PFEGammaAlgo::_refinableObjects
std::list< ProtoEGObject > _refinableObjects
Definition: PFEGammaAlgo.h:201
first
bool first
Definition: L1TdeRCT.cc:75
reco::PFBlockElementGsfTrack::trackType
virtual bool trackType(TrackType trType) const
Definition: PFBlockElementGsfTrack.h:33
cond::ecalcond::value_type
Container::value_type value_type
Definition: EcalChannelStatusPyWrapper.cc:33
PFEGammaAlgo::outcandsextra_
reco::PFCandidateEGammaExtraCollection outcandsextra_
Definition: PFEGammaAlgo.h:177
PFEGammaAlgo::linkKFTrackToECAL
void linkKFTrackToECAL(const PFKFFlaggedElement &, ProtoEGObject &)
Definition: PFEGammaAlgo.cc:1683
PFEGammaAlgo::buildAndRefineEGObjects
void buildAndRefineEGObjects(const reco::PFBlockRef &block)
Definition: PFEGammaAlgo.cc:849
reco::PFCandidate::setEcalEnergy
void setEcalEnergy(float eeRaw, float eeCorr)
set corrected Ecal energy
Definition: PFCandidate.h:200
PFLayer::Layer
Layer
layer definition
Definition: PFLayer.h:31
PFEGammaAlgo::_currentlinks
reco::PFBlock::LinkData _currentlinks
Definition: PFEGammaAlgo.h:184
reco::e1
Float e1
Definition: deltaR.h:22
reco::PFCandidateCollection
std::vector< reco::PFCandidate > PFCandidateCollection
collection of PFCandidates
Definition: PFCandidateFwd.h:12
reco::Vertex::x
double x() const
x coordinate
Definition: Vertex.h:108
edm::Ptr::key
key_type key() const
Definition: Ptr.h:169
reco::PFBlockElementTrack::trackType
virtual bool trackType(TrackType trType) const
Definition: PFBlockElementTrack.h:30
reco::SuperCluster::rawEnergy
double rawEnergy() const
raw uncorrected energy (sum of energies of component BasicClusters)
Definition: SuperCluster.h:47
PFEGammaAlgo::fill_extra_info
void fill_extra_info(const ProtoEGObject &, reco::PFCandidateEGammaExtra &)
Definition: PFEGammaAlgo.cc:2195
PFEGammaAlgo::egCandidate_
reco::PFCandidateCollection egCandidate_
Definition: PFEGammaAlgo.h:380
patCandidatesForDimuonsSequences_cff.ecal
tuple ecal
Definition: patCandidatesForDimuonsSequences_cff.py:70
SiPixelRawToDigiRegional_cfi.deltaPhi
tuple deltaPhi
Definition: SiPixelRawToDigiRegional_cfi.py:9
PFEGammaAlgo::isAMuon
bool isAMuon(const reco::PFBlockElement &)
Definition: PFEGammaAlgo.cc:822
reco::PFCandidateEGammaExtra::setDeltaEta
void setDeltaEta(float val)
set the delta eta
Definition: PFCandidateEGammaExtra.cc:90
reco::PFCandidate::setGsfTrackRef
void setGsfTrackRef(const reco::GsfTrackRef &ref)
set gsftrack reference
Definition: PFCandidate.cc:450
PFEGammaAlgo::ProtoEGObject::electronSeed
reco::ElectronSeedRef electronSeed
Definition: PFEGammaAlgo.h:79
PFEGammaAlgo::calculate_ele_mva
float calculate_ele_mva(const ProtoEGObject &, reco::PFCandidateEGammaExtra &)
Definition: PFEGammaAlgo.cc:2031
math::XYZVector
XYZVectorD XYZVector
spatial vector with cartesian internal representation
Definition: Vector3D.h:30
reco::PFBlock::associatedElements
void associatedElements(unsigned i, const LinkData &linkData, std::multimap< double, unsigned > &sortedAssociates, reco::PFBlockElement::Type type=PFBlockElement::NONE, LinkTest test=LINKTEST_RECHIT) const
Definition: PFBlock.cc:75
math::XYZPoint
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
customizeTrackingMonitorSeedNumber.idx
tuple idx
DEBUGGING if hasattr(process,&quot;trackMonIterativeTracking2012&quot;): print &quot;trackMonIterativeTracking2012 D...
Definition: customizeTrackingMonitorSeedNumber.py:15
reco::SuperCluster::addPreshowerCluster
void addPreshowerCluster(const CaloClusterPtr &r)
add reference to constituent BasicCluster
Definition: SuperCluster.h:117
PFEGammaAlgo::SigmaEtaEta
float SigmaEtaEta
Definition: PFEGammaAlgo.h:311
run_regression.test
tuple test
Definition: run_regression.py:386
b
double b
Definition: hdecay.h:120
docast
#define docast(x, y)
Definition: PFEGammaAlgo.cc:42
PFEGammaAlgo::linkRefinableObjectGSFTracksToKFs
void linkRefinableObjectGSFTracksToKFs(ProtoEGObject &)
Definition: PFEGammaAlgo.cc:1523
PFEGammaAlgo::convGsfTrack_
std::vector< std::pair< unsigned int, unsigned int > > convGsfTrack_
Definition: PFEGammaAlgo.h:290
reco::PFCandidateEGammaExtra::setSuperClusterPFECALRef
void setSuperClusterPFECALRef(reco::SuperClusterRef sc)
set reference to the corresponding supercluster
Definition: PFCandidateEGammaExtra.h:114
reco::PFCandidateEGammaExtra::addConversionRef
void addConversionRef(const reco::ConversionRef &convref)
add Conversions from PF
Definition: PFCandidateEGammaExtra.cc:160
PFEGammaAlgo::tmvaReaderEle_
TMVA::Reader * tmvaReaderEle_
Definition: PFEGammaAlgo.h:294
PFEGammaAlgo::ProtoEGObject::primaryKFs
std::vector< PFKFFlaggedElement > primaryKFs
Definition: PFEGammaAlgo.h:88
reco::e2
Float e2
Definition: deltaR.h:23
PFEGammaAlgo::ClusterMap
std::unordered_map< const PFClusterElement *, std::vector< PFClusterFlaggedElement > > ClusterMap
Definition: PFEGammaAlgo.h:71
GetRecoTauVFromDQM_MC_cff.cl
tuple cl
Definition: GetRecoTauVFromDQM_MC_cff.py:37
PFEGammaAlgo::EtotPinMode
float EtotPinMode
Definition: PFEGammaAlgo.h:309
reco::SuperCluster::addCluster
void addCluster(const CaloClusterPtr &r)
add reference to constituent BasicCluster
Definition: SuperCluster.h:111
reco::PFBlockElementCluster::Dump
void Dump(std::ostream &out=std::cout, const char *tab=" ") const
print the object inside the element
Definition: PFBlockElementCluster.cc:11
PFEGammaAlgo::_recoveredlinks
ElementMap _recoveredlinks
Definition: PFEGammaAlgo.h:188
PFEGammaAlgo::removeOrLinkECALClustersToKFTracks
void removeOrLinkECALClustersToKFTracks()
Definition: PFEGammaAlgo.cc:1358
reco::PFCandidate
Particle reconstructed by the particle flow algorithm.
Definition: PFCandidate.h:38
begin
#define begin
Definition: vmac.h:30
reco::PFBlockElementSuperCluster::Dump
void Dump(std::ostream &out=std::cout, const char *tab=" ") const
print the object inside the element
Definition: PFBlockElementSuperCluster.cc:11
combine.key
list key
Definition: combine.py:13
LOGWARN
#define LOGWARN(x)
Definition: PFEGammaAlgo.cc:44
PFEGammaAlgo::buildRefinedSuperCluster
reco::SuperCluster buildRefinedSuperCluster(const ProtoEGObject &)
Definition: PFEGammaAlgo.cc:2222
PFEGammaAlgo::X0_inner
TH2D * X0_inner
Definition: PFEGammaAlgo.h:371
PFEGammaAlgo::unlinkRefinableObjectKFandECALMatchedToHCAL
void unlinkRefinableObjectKFandECALMatchedToHCAL(ProtoEGObject &, bool removeFreeECAL=false, bool removeSCECAL=false)
Definition: PFEGammaAlgo.cc:2425
reco::LeafCandidate::setP4
virtual void setP4(const LorentzVector &p4) GCC11_FINAL
set 4-momentum
Definition: LeafCandidate.h:190
a
double a
Definition: hdecay.h:121
PFEGammaAlgo::PFKFFlaggedElement
std::pair< const PFKFElement *, bool > PFKFFlaggedElement
Definition: PFEGammaAlgo.h:63
reco::PFBlockElementTrack::trackRefPF
const PFRecTrackRef & trackRefPF() const
Definition: PFBlockElementTrack.h:54
PFEGammaAlgo::mergeROsByAnyLink
void mergeROsByAnyLink(std::list< ProtoEGObject > &)
Definition: PFEGammaAlgo.cc:1440
reco::PFCandidate::setSuperClusterRef
void setSuperClusterRef(const reco::SuperClusterRef &scRef)
Definition: PFCandidate.cc:616
PFEGammaAlgo::attachPSClusters
int attachPSClusters(const PFClusterElement *, ClusterMap::mapped_type &)
Definition: PFEGammaAlgo.cc:1253
relval_parameters_module.energy
string energy
Definition: relval_parameters_module.py:29
edm::Ref::index
key_type index() const
Definition: Ref.h:269
PFEGammaAlgo::linkRefinableObjectKFTracksToECAL
void linkRefinableObjectKFTracksToECAL(ProtoEGObject &)
Definition: PFEGammaAlgo.cc:1676
reco::PFCandidateEGammaExtra::setKfTrackRef
void setKfTrackRef(const reco::TrackRef &ref)
set kf track reference
Definition: PFCandidateEGammaExtra.cc:56
reco::SuperCluster::seed
const CaloClusterPtr & seed() const
seed BasicCluster
Definition: SuperCluster.h:66
reco::PFCandidateEGammaExtra::setMVA
void setMVA(float val)
set the result (mostly for debugging)
Definition: PFCandidateEGammaExtra.cc:124
edm::false
volatile std::atomic< bool > shutdown_flag false
Definition: UnixSignalHandlers.cc:22
PFEGammaAlgo::tmvaReader_
TMVA::Reader * tmvaReader_
Definition: PFEGammaAlgo.h:330
PFEGammaAlgo::linkRefinableObjectSecondaryKFsToECAL
void linkRefinableObjectSecondaryKFsToECAL(ProtoEGObject &)
Definition: PFEGammaAlgo.cc:1881
reco::PFCandidateEGammaExtra::setLateBrem
void setLateBrem(float val)
set LateBrem
Definition: PFCandidateEGammaExtra.cc:69
reco::Vertex::tracks_begin
trackRef_iterator tracks_begin() const
first iterator over tracks
Definition: Vertex.cc:39
PFEGammaAlgo::_finalCandidates
reco::PFCandidateCollection _finalCandidates
Definition: PFEGammaAlgo.h:203
PFEGammaAlgo::X0_middle
TH2D * X0_middle
Definition: PFEGammaAlgo.h:372
PFEGammaAlgo::linkRefinableObjectConvSecondaryKFsToSecondaryKFs
void linkRefinableObjectConvSecondaryKFsToSecondaryKFs(ProtoEGObject &)
Definition: PFEGammaAlgo.cc:1814
PFEGammaAlgo::verbosityLevel_
verbosityLevel verbosityLevel_
Definition: PFEGammaAlgo.h:323
edm::Ref::get
T const * get() const
Returns C++ pointer to the item.
Definition: Ref.h:242
reco::PFCandidate::setTrackRef
void setTrackRef(const reco::TrackRef &ref)
set track reference
Definition: PFCandidate.cc:412
reco::LeafCandidate::setPdgId
virtual void setPdgId(int pdgId) GCC11_FINAL
Definition: LeafCandidate.h:238
PFEGammaAlgo::outcands_
reco::PFCandidateCollection outcands_
Definition: PFEGammaAlgo.h:176
PFEGammaAlgo::fillPFCandidates
void fillPFCandidates(const std::list< ProtoEGObject > &, reco::PFCandidateCollection &, reco::PFCandidateEGammaExtraCollection &)
Definition: PFEGammaAlgo.cc:1904
stdcomb::next_combination
bool next_combination(BidIt n_begin, BidIt n_end, BidIt r_begin, BidIt r_end)
Definition: combination.h:22
patCandidatesForDimuonsSequences_cff.hcal
tuple hcal
Definition: patCandidatesForDimuonsSequences_cff.py:71
reco::SuperCluster::setPreshowerEnergyPlane1
void setPreshowerEnergyPlane1(double preshowerEnergy1)
Definition: SuperCluster.h:60
PFEGammaAlgo::ProtoEGObject::primaryGSFs
std::vector< PFGSFFlaggedElement > primaryGSFs
Definition: PFEGammaAlgo.h:86
findQualityFiles.size
tuple size
Write out results.
Definition: findQualityFiles.py:442
Algo
Definition: fakeMenu.h:6
constexpr
#define constexpr
Definition: GCC11Compatibility.h:36
reco::SuperCluster::preshowerClustersEnd
CaloCluster_iterator preshowerClustersEnd() const
last iterator over PreshowerCluster constituents
Definition: SuperCluster.h:84
reco::PFCandidateEGammaExtra::addExtraNonConvTrack
void addExtraNonConvTrack(const reco::PFBlockRef &blk, const reco::PFBlockElementTrack &tkref)
track counting for electrons and photons
Definition: PFCandidateEGammaExtra.h:180
reco::SuperCluster::setPreshowerEnergy
void setPreshowerEnergy(double preshowerEnergy)
Definition: SuperCluster.h:59
PFEGammaAlgo::linkRefinableObjectPrimaryGSFTrackToHCAL
void linkRefinableObjectPrimaryGSFTrackToHCAL(ProtoEGObject &)
Definition: PFEGammaAlgo.cc:1650
reco::PFBlock
Block of elements.
Definition: PFBlock.h:30
Generated for CMSSW Reference Manual by doxygen 1.8.5