CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
List of all members | Public Member Functions | Private Attributes
RBXAndHPDCleaner Class Reference

#include <RBXAndHPDCleaner.h>

Inheritance diagram for RBXAndHPDCleaner:
RecHitTopologicalCleanerBase

Public Member Functions

void clean (const edm::Handle< reco::PFRecHitCollection > &input, std::vector< bool > &mask)
 
RBXAndHPDCleaner & operator= (const RBXAndHPDCleaner &)=delete
 
 RBXAndHPDCleaner (const edm::ParameterSet &conf)
 
 RBXAndHPDCleaner (const RBXAndHPDCleaner &)=delete
 
- Public Member Functions inherited from RecHitTopologicalCleanerBase
const std::string & name () const
 
RecHitTopologicalCleanerBase & operator= (const RecHitTopologicalCleanerBase &)=delete
 
 RecHitTopologicalCleanerBase (const edm::ParameterSet &conf)
 
 RecHitTopologicalCleanerBase (const RecHitTopologicalCleanerBase &)=delete
 

Private Attributes

std::unordered_map< int,
std::vector< unsigned > > 		_hpds
 
std::unordered_map< int,
std::vector< unsigned > > 		_rbxs
 

Detailed Description

Definition at line 8 of file RBXAndHPDCleaner.h.

Constructor & Destructor Documentation

RBXAndHPDCleaner::RBXAndHPDCleaner ( const edm::ParameterSet &  conf)
inline

Definition at line 10 of file RBXAndHPDCleaner.h.

10  :
11  RecHitTopologicalCleanerBase(conf) { }
RecHitTopologicalCleanerBase::RecHitTopologicalCleanerBase
RecHitTopologicalCleanerBase(const edm::ParameterSet &conf)
Definition: RecHitTopologicalCleanerBase.h:13
RBXAndHPDCleaner::RBXAndHPDCleaner ( const RBXAndHPDCleaner &  )
delete

Member Function Documentation

void RBXAndHPDCleaner::clean ( const edm::Handle< reco::PFRecHitCollection > &  input,
std::vector< bool > &  mask 
)
virtual

Implements RecHitTopologicalCleanerBase.

Definition at line 15 of file RBXAndHPDCleaner.cc.

References _hpds, _rbxs, funct::abs(), reco::PFRecHit::detId(), alignCSCRings::e, reco::PFRecHit::energy(), Basic3DVector< T >::eta(), PFLayer::HCAL_BARREL1, PFLayer::HCAL_ENDCAP, i, customizeTrackingMonitorSeedNumber::idx, HcalDetId::ieta(), HcalDetId::iphi(), relval_2017::k, reco::PFRecHit::layer(), reco::PFRecHit::neighbours4(), Basic3DVector< T >::phi(), reco::PFRecHit::position(), HI_PhotonSkim_cff::rechits, mathSSE::sqrt(), and dtDQMClient_cfg::threshold.

16  {
17  _hpds.clear(); _rbxs.clear();
18  //need to run over energy sorted rechits
19  std::vector<std::pair<unsigned,double> > ordered_hits;
20  ordered_hits.reserve(input->size());
21  for( unsigned i = 0; i < input->size(); ++i ) {
22  std::pair<unsigned,double> val = std::make_pair(i,input->at(i).energy());
23  auto pos = std::upper_bound(ordered_hits.begin(),ordered_hits.end(),
24  val, greaterByEnergy);
25  ordered_hits.insert(pos,val);
26  }
27 
28  for( const auto& idx_e : ordered_hits ) {
29  if( !mask[idx_e.first] ) continue;
30  const unsigned idx = idx_e.first;
31  const reco::PFRecHit& rechit = input->at(idx);
32  int layer = rechit.layer();
33  if ( layer != PFLayer::HCAL_BARREL1 && layer != PFLayer::HCAL_ENDCAP ) {
34  break;
35  }
36  HcalDetId theHcalDetId(rechit.detId());
37  int ieta = theHcalDetId.ieta();
38  int iphi = theHcalDetId.iphi();
39  int ihpd = 0, irbx = 0;
40  switch( layer ) {
41  case PFLayer::HCAL_BARREL1:
42  ihpd = ( ieta < 0 ? -iphi : iphi );
43  irbx = ( ieta < 0 ? -(iphi+5)/4 : (iphi+5)/4 );
44  break;
45  case PFLayer::HCAL_ENDCAP:
46  ihpd = ( ieta < 0 ? -(iphi+1)/2-100 : (iphi+1)/2+100 );
47  irbx = ( ieta < 0 ? -(iphi+5)/4-20 : (iphi+5)/4+20 );
48  break;
49  default:
50  break;
51  }
52  switch( std::abs(irbx) ) {
53  case 19:
54  irbx = ( irbx < 0 ? -1 : 1 );
55  break;
56  case 39:
57  irbx = ( irbx < 0 ? -21 : 21 );
58  break;
59  default:
60  break;
61  }
62  _hpds[ihpd].push_back(idx);
63  _rbxs[irbx].push_back(idx);
64  }
65  // loop on the rbx's we found and clean RBX's with tons of rechits
66  // and lots of energy
67  double totalEta = 0., totalEtaW = 0., totalPhi = 0., totalPhiW = 0.,
68  totalEnergy = 0.;
69  double totalEta2 = 1E-9, totalEta2W = 1E-9, totalPhi2 = 1E-9,
70  totalPhi2W = 1E-9, totalEnergy2 = 1E-9;
71  unsigned nSeeds = 0, nSeeds0 = 0;
72  std::unordered_map<int, std::vector<unsigned> > theHPDs;
73  std::unordered_multimap<double, unsigned> theEnergies;
74  for( const auto& itrbx : _rbxs ) {
75  if( ( std::abs(itrbx.first)<20 && itrbx.second.size() > 30 ) ||
76  ( std::abs(itrbx.first)>20 && itrbx.second.size() > 30 ) ) {
77  const std::vector<unsigned>& rechits = itrbx.second;
78  theHPDs.clear();
79  theEnergies.clear();
80  totalEta = totalEtaW = totalPhi = totalPhiW = totalEnergy = 0.;
81  totalEta2 = totalEta2W = totalPhi2 = totalPhi2W = totalEnergy2 = 1e-9;
82  nSeeds = nSeeds0 = rechits.size();
83  for( unsigned jh = 0; jh < rechits.size(); ++jh ) {
84  const reco::PFRecHit& rechit = (*input)[jh];
85  // check if rechit is a seed
86  unsigned nN = 0 ; // neighbours over threshold
87  bool isASeed = true;
88  auto const & neighbours4 = rechit.neighbours4();
89  for( auto k : neighbours4 ) {
90  auto const & neighbour = (*input)[k];
91  if( neighbour.energy() > rechit.energy() ) {
92  --nSeeds; --nSeeds0;
93  isASeed = false;
94  break;
95  } else {
96  if( neighbour.energy() > 0.4 ) ++nN;
97  }
98  }
99  if ( isASeed && !nN ) --nSeeds0;
100 
101  HcalDetId theHcalDetId(rechit.detId());
102  int iphi = theHcalDetId.iphi();
103  switch( rechit.layer() ) {
104  case PFLayer::HCAL_BARREL1:
105  theHPDs[iphi].push_back(rechits[jh]);
106  break;
107  case PFLayer::HCAL_ENDCAP:
108  theHPDs[(iphi-1)/2].push_back(rechits[jh]);
109  break;
110  default:
111  break;
112  }
113  const double rhenergy = rechit.energy();
114  const double rhphi = rechit.position().phi();
115  const double rhphi2 = rhphi*rhphi;
116  const double rheta = rechit.position().eta();
117  const double rheta2 = rheta*rheta;
118  theEnergies.emplace(rhenergy,rechits[jh]);
119  totalEnergy += rhenergy;
120  totalPhi += std::abs(rhphi);
121  totalPhiW += std::abs(rhphi)*rhenergy;
122  totalEta += rheta;
123  totalEtaW += rheta*rhenergy;
124  totalEnergy2 += rhenergy*rhenergy;
125  totalPhi2 += rhphi2;
126  totalPhi2W += rhphi2*rhenergy;
127  totalEta2 += rheta2;
128  totalEta2W += rheta2*rhenergy;
129  }
130  totalPhi /= rechits.size();
131  totalEta /= rechits.size();
132  totalPhiW /= totalEnergy;
133  totalEtaW /= totalEnergy;
134  totalPhi2 /= rechits.size();
135  totalEta2 /= rechits.size();
136  totalPhi2W /= totalEnergy;
137  totalEta2W /= totalEnergy;
138  totalPhi2 = std::sqrt(totalPhi2 - totalPhi*totalPhi);
139  totalEta2 = std::sqrt(totalEta2 - totalEta*totalEta);
140  totalPhi2W = std::sqrt(totalPhi2W - totalPhiW*totalPhi2);
141  totalEta2W = std::sqrt(totalEta2W - totalEtaW*totalEtaW);
142  totalEnergy /= rechits.size();
143  totalEnergy2 /= rechits.size();
144  totalEnergy2 = std::sqrt(totalEnergy2 - totalEnergy*totalEnergy);
145  if( nSeeds0 > 6 ) {
146  unsigned nHPD15 = 0;
147  for( const auto& itHPD : theHPDs ) {
148  int hpdN = itHPD.first;
149  const std::vector<unsigned>& hpdHits = itHPD.second;
150  if( ( std::abs(hpdN) < 100 && hpdHits.size() > 14 ) ||
151  ( std::abs(hpdN) > 100 && hpdHits.size() > 14 ) ) ++nHPD15;
152  }
153  if( nHPD15 > 1 ) {
154  unsigned nn = 0;
155  double threshold = 1.0;
156  for( const auto& itEn : theEnergies ) {
157  ++nn;
158  if( nn< 5 ) {
159  mask[itEn.second] = false;
160  } else if ( nn == 5 ) {
161  threshold = itEn.first*5;
162  mask[itEn.second] = false;
163  } else {
164  if( itEn.first < threshold ) mask[itEn.second] = false;
165  }
166  }
167  }
168  }
169  }
170  }
171 
172  // loop on the HPDs
173  std::unordered_map<int,std::vector<unsigned> >::iterator neighbour1;
174  std::unordered_map<int,std::vector<unsigned> >::iterator neighbour2;
175  std::unordered_map<int,std::vector<unsigned> >::iterator neighbour0;
176  std::unordered_map<int,std::vector<unsigned> >::iterator neighbour3;
177  unsigned size1 = 0, size2 = 0;
178  for( const auto& ithpd : _hpds ) {
179  const std::vector<unsigned>& rechits = ithpd.second;
180  theEnergies.clear();
181  totalEnergy = 0;
182  totalEnergy2 = 1e-9;
183  for( const unsigned rhidx : rechits ) {
184  const reco::PFRecHit & rechit = input->at(rhidx);
185  const double e = rechit.energy();
186  totalEnergy += e;
187  totalEnergy2 += e*e;
188  theEnergies.emplace(rechit.energy(),rhidx);
189  }
190  totalEnergy /= rechits.size();
191  totalEnergy2 /= rechits.size();
192  totalEnergy2 = std::sqrt(totalEnergy2 - totalEnergy*totalEnergy);
193 
194  const int thehpd = ithpd.first;
195  switch( std::abs(thehpd) ) {
196  case 1:
197  neighbour1 = ( thehpd > 0 ? _hpds.find(72) : _hpds.find(-72) );
198  break;
199  case 72:
200  neighbour2 = ( thehpd > 0 ? _hpds.find(1) : _hpds.find(-1) );
201  break;
202  case 101:
203  neighbour1 = ( thehpd > 0 ? _hpds.find(136) : _hpds.find(-136) );
204  break;
205  case 136:
206  neighbour2 = ( thehpd > 0 ? _hpds.find(101) : _hpds.find(-101) );
207  break;
208  default:
209  neighbour1 = ( thehpd > 0 ? _hpds.find(thehpd-1) : _hpds.find(thehpd+1) );
210  neighbour2 = ( thehpd > 0 ? _hpds.find(thehpd+1) : _hpds.find(thehpd-1) );
211  break;
212  }
213  if( neighbour1 != _hpds.end() ) {
214  const int nb1 = neighbour1->first;
215  switch( std::abs(nb1) ) {
216  case 1:
217  neighbour0 = ( nb1 > 0 ? _hpds.find(72) : _hpds.find(-72) );
218  break;
219  case 101:
220  neighbour0 = ( nb1 > 0 ? _hpds.find(136) : _hpds.find(-136) );
221  break;
222  default:
223  neighbour0 = ( nb1 > 0 ? _hpds.find(nb1-1) : _hpds.find(nb1+1) );
224  break;
225  }
226  } else {
227  neighbour0 = _hpds.end();
228  }
229 
230  if( neighbour2 != _hpds.end() ) {
231  const int nb2 = neighbour2->first;
232  switch( std::abs(nb2) ) {
233  case 72:
234  neighbour3 = ( nb2 > 0 ? _hpds.find(1) : _hpds.find(-1) );
235  break;
236  case 136:
237  neighbour3 = ( nb2 > 0 ? _hpds.find(101) : _hpds.find(-101) );
238  break;
239  default:
240  neighbour3 = ( nb2 > 0 ? _hpds.find(nb2+1) : _hpds.find(nb2-1) );
241  break;
242  }
243  } else {
244  neighbour3 = _hpds.end();
245  }
246 
247  size1 = neighbour1 != _hpds.end() ? neighbour1->second.size() : 0;
248  size2 = neighbour2 != _hpds.end() ? neighbour2->second.size() : 0;
249  if( size1 > 10 ) {
250  if ( ( abs(neighbour1->first) > 100 && neighbour1->second.size() > 15 ) ||
251  ( abs(neighbour1->first) < 100 && neighbour1->second.size() > 12 ) )
252  size1 = neighbour0 != _hpds.end() ? neighbour0->second.size() : 0;
253  }
254  if( size2 > 10 ) {
255  if ( ( abs(neighbour2->first) > 100 && neighbour2->second.size() > 15 ) ||
256  ( abs(neighbour2->first) < 100 && neighbour2->second.size() > 12 ) )
257  size2 = neighbour3 != _hpds.end() ? neighbour3->second.size() : 0;
258  }
259  if( ( std::abs(ithpd.first) > 100 && ithpd.second.size() > 15 ) ||
260  ( std::abs(ithpd.first) < 100 && ithpd.second.size() > 12 ) ) {
261  if( (double)(size1+size2)/(float)ithpd.second.size() < 1.0 ) {
262  unsigned nn = 0;
263  double threshold = 1.0;
264  for( const auto& itEn : theEnergies ) {
265  if( nn < 5 ) {
266  mask[itEn.second] = false;
267  } else if ( nn == 5 ) {
268  threshold = itEn.first*2.5;
269  mask[itEn.second] = false;
270  } else {
271  if( itEn.first < threshold ) mask[itEn.second] = false;
272  }
273  }
274  }
275  }
276  }
277 }
i
int i
Definition: DBlmapReader.cc:9
reco::PFRecHit::detId
unsigned detId() const
rechit detId
Definition: PFRecHit.h:108
relval_2017.k
list k
Definition: relval_2017.py:23
Basic3DVector::phi
Geom::Phi< T > phi() const
Definition: extBasic3DVector.h:137
reco::PFRecHit::position
PositionType const & position() const
rechit cell centre x, y, z
Definition: PFRecHit.h:129
reco::PFRecHit::layer
PFLayer::Layer layer() const
rechit layer
Definition: PFRecHit.h:111
reco::PFRecHit
Particle flow rechit (rechit + geometry and topology information). See clustering algorithm in PFClus...
Definition: PFRecHit.h:31
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
HcalDetId::ieta
int ieta() const
get the cell ieta
Definition: HcalDetId.h:56
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
dtDQMClient_cfg.threshold
tuple threshold
Definition: dtDQMClient_cfg.py:16
reco::PFRecHit::energy
float energy() const
rechit energy
Definition: PFRecHit.h:114
RBXAndHPDCleaner::_hpds
std::unordered_map< int, std::vector< unsigned > > _hpds
Definition: RBXAndHPDCleaner.h:19
RBXAndHPDCleaner::_rbxs
std::unordered_map< int, std::vector< unsigned > > _rbxs
Definition: RBXAndHPDCleaner.h:19
HcalDetId::iphi
int iphi() const
get the cell iphi
Definition: HcalDetId.cc:124
customizeTrackingMonitorSeedNumber.idx
tuple idx
DEBUGGING if hasattr(process,&quot;trackMonIterativeTracking2012&quot;): print &quot;trackMonIterativeTracking2012 D...
Definition: customizeTrackingMonitorSeedNumber.py:15
HI_PhotonSkim_cff.rechits
tuple rechits
Definition: HI_PhotonSkim_cff.py:76
reco::PFRecHit::neighbours4
Neighbours neighbours4() const
Definition: PFRecHit.h:87
RBXAndHPDCleaner& RBXAndHPDCleaner::operator= ( const RBXAndHPDCleaner &  )
delete

Member Data Documentation

std::unordered_map<int,std::vector<unsigned> > RBXAndHPDCleaner::_hpds
private

Definition at line 19 of file RBXAndHPDCleaner.h.

Referenced by clean().

std::unordered_map<int,std::vector<unsigned> > RBXAndHPDCleaner::_rbxs
private

Definition at line 19 of file RBXAndHPDCleaner.h.

Referenced by clean().

Generated for CMSSW Reference Manual by doxygen 1.8.5