CMS 3D CMS Logo

LocalMaximumSeedFinder.cc
Go to the documentation of this file.
2 
3 #include <algorithm>
4 #include <queue>
5 #include <cfloat>
8 
9 namespace {
10  const reco::PFRecHit::Neighbours _noNeighbours(nullptr,0);
11 }
12 
15  SeedFinderBase(conf),
16  _nNeighbours(conf.getParameter<int>("nNeighbours")),
17  _layerMap({ {"PS2",(int)PFLayer::PS2},
18  {"PS1",(int)PFLayer::PS1},
19  {"ECAL_ENDCAP",(int)PFLayer::ECAL_ENDCAP},
20  {"ECAL_BARREL",(int)PFLayer::ECAL_BARREL},
21  {"NONE",(int)PFLayer::NONE},
22  {"HCAL_BARREL1",(int)PFLayer::HCAL_BARREL1},
23  {"HCAL_BARREL2_RING0",(int)PFLayer::HCAL_BARREL2},
24  // hack to deal with ring1 in HO
25  {"HCAL_BARREL2_RING1", 19},
26  {"HCAL_ENDCAP",(int)PFLayer::HCAL_ENDCAP},
27  {"HF_EM",(int)PFLayer::HF_EM},
28  {"HF_HAD",(int)PFLayer::HF_HAD} }) {
29  const std::vector<edm::ParameterSet>& thresholds =
30  conf.getParameterSetVector("thresholdsByDetector");
31  for( const auto& pset : thresholds ) {
32  const std::string& det = pset.getParameter<std::string>("detector");
33  const double& thresh_E = pset.getParameter<double>("seedingThreshold");
34  const double& thresh_pT = pset.getParameter<double>("seedingThresholdPt");
35  const double thresh_pT2 = thresh_pT*thresh_pT;
36  auto entry = _layerMap.find(det);
37  if( entry == _layerMap.end() ) {
38  throw cms::Exception("InvalidDetectorLayer")
39  << "Detector layer : " << det << " is not in the list of recognized"
40  << " detector layers!";
41  }
42  _thresholds[entry->second+layerOffset]=
43  std::make_pair(thresh_E,thresh_pT2);
44  }
45 }
46 
47 // the starting state of seedable is all false!
50  const std::vector<bool>& mask,
51  std::vector<bool>& seedable ) {
52 
53  auto nhits = input->size();
54  initDynArray(bool,nhits,usable,true);
55  //need to run over energy sorted rechits
56  declareDynArray(float,nhits,energies);
57  unInitDynArray(int,nhits,qst); // queue storage
58  auto cmp = [&](int i, int j) { return energies[i] < energies[j]; };
59  std::priority_queue<int, DynArray<int>, decltype(cmp)> ordered_hits(cmp,std::move(qst));
60 
61  for( unsigned i = 0; i < nhits; ++i ) {
62  if( !mask[i] ) continue; // cannot seed masked objects
63  const auto & maybeseed = (*input)[i];
64  energies[i]=maybeseed.energy();
65  int seedlayer = (int)maybeseed.layer();
66  if( seedlayer == PFLayer::HCAL_BARREL2 &&
67  std::abs(maybeseed.positionREP().eta()) > 0.34 ) {
68  seedlayer = 19;
69  }
70  auto const & thresholds = _thresholds[seedlayer+layerOffset];
71  if( maybeseed.energy() < thresholds.first ||
72  maybeseed.pt2() < thresholds.second ) usable[i] = false;
73  if( !usable[i] ) continue;
74  ordered_hits.push(i);
75  }
76 
77 
78  while(!ordered_hits.empty() ) {
79  auto idx = ordered_hits.top();
80  ordered_hits.pop();
81  if( !usable[idx] ) continue;
82  //get the neighbours of this seed
83  const auto & maybeseed = (*input)[idx];
84  reco::PFRecHit::Neighbours myNeighbours;
85  switch( _nNeighbours ) {
86  case -1:
87  myNeighbours = maybeseed.neighbours();
88  break;
89  case 0: // for HF clustering
90  myNeighbours = _noNeighbours;
91  break;
92  case 4:
93  myNeighbours = maybeseed.neighbours4();
94  break;
95  case 8:
96  myNeighbours = maybeseed.neighbours8();
97  break;
98  default:
99  throw cms::Exception("InvalidConfiguration")
100  << "LocalMaximumSeedFinder only accepts nNeighbors = {-1,0,4,8}";
101  }
102  seedable[idx] = true;
103  for( auto neighbour : myNeighbours ) {
104  if( !mask[neighbour] ) continue;
105  if( energies[neighbour] > energies[idx] ) {
106 // std::cout << "how this can be?" << std::endl;
107  seedable[idx] = false;
108  break;
109  }
110  }
111  if( seedable[idx] ) {
112  for( auto neighbour : myNeighbours ) {
113  usable[neighbour] = false;
114  }
115  }
116  }
117 
118  LogDebug("LocalMaximumSeedFinder") << " found " << std::count(seedable.begin(),seedable.end(),true) << " seeds";
119 
120 }
#define LogDebug(id)
const std::unordered_map< std::string, int > _layerMap
#define initDynArray(T, n, x, i)
Definition: DynArray.h:60
void findSeeds(const edm::Handle< reco::PFRecHitCollection > &input, const std::vector< bool > &mask, std::vector< bool > &seedable)
LocalMaximumSeedFinder(const edm::ParameterSet &conf)
static std::string const input
Definition: EdmProvDump.cc:44
std::array< std::pair< double, double >, 35 > _thresholds
#define unInitDynArray(T, n, x)
Definition: DynArray.h:58
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
#define declareDynArray(T, n, x)
Definition: DynArray.h:59
def move(src, dest)
Definition: eostools.py:510