test
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(), PFLayer::HCAL_BARREL1, PFLayer::HCAL_ENDCAP, i, customizeTrackingMonitorSeedNumber::idx, HcalDetId::ieta(), HcalDetId::iphi(), reco::PFRecHit::layer(), reco::PFRecHit::neighbours4(), 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->at(jh);
85  // check if rechit is a seed
86  unsigned nN = 0 ; // neighbours over threshold
87  bool isASeed = true;
88  const reco::PFRecHitRefVector& neighbours4 = rechit.neighbours4();
89  for( const reco::PFRecHitRef& neighbour : neighbours4 ) {
90  if( neighbour->energy() > rechit.energy() ) {
91  --nSeeds; --nSeeds0;
92  isASeed = false;
93  break;
94  } else {
95  if( neighbour->energy() > 0.4 ) ++nN;
96  }
97  }
98  if ( isASeed && !nN ) --nSeeds0;
99 
100  HcalDetId theHcalDetId(rechit.detId());
101  int iphi = theHcalDetId.iphi();
102  switch( rechit.layer() ) {
103  case PFLayer::HCAL_BARREL1:
104  theHPDs[iphi].push_back(rechits[jh]);
105  break;
106  case PFLayer::HCAL_ENDCAP:
107  theHPDs[(iphi-1)/2].push_back(rechits[jh]);
108  break;
109  default:
110  break;
111  }
112  const double rhenergy = rechit.energy();
113  const double rhphi = rechit.position().phi();
114  const double rhphi2 = rhphi*rhphi;
115  const double rheta = rechit.position().eta();
116  const double rheta2 = rheta*rheta;
117  theEnergies.emplace(rhenergy,rechits[jh]);
118  totalEnergy += rhenergy;
119  totalPhi += std::abs(rhphi);
120  totalPhiW += std::abs(rhphi)*rhenergy;
121  totalEta += rheta;
122  totalEtaW += rheta*rhenergy;
123  totalEnergy2 += rhenergy*rhenergy;
124  totalPhi2 += rhphi2;
125  totalPhi2W += rhphi2*rhenergy;
126  totalEta2 += rheta2;
127  totalEta2W += rheta2*rhenergy;
128  }
129  totalPhi /= rechits.size();
130  totalEta /= rechits.size();
131  totalPhiW /= totalEnergy;
132  totalEtaW /= totalEnergy;
133  totalPhi2 /= rechits.size();
134  totalEta2 /= rechits.size();
135  totalPhi2W /= totalEnergy;
136  totalEta2W /= totalEnergy;
137  totalPhi2 = std::sqrt(totalPhi2 - totalPhi*totalPhi);
138  totalEta2 = std::sqrt(totalEta2 - totalEta*totalEta);
139  totalPhi2W = std::sqrt(totalPhi2W - totalPhiW*totalPhi2);
140  totalEta2W = std::sqrt(totalEta2W - totalEtaW*totalEtaW);
141  totalEnergy /= rechits.size();
142  totalEnergy2 /= rechits.size();
143  totalEnergy2 = std::sqrt(totalEnergy2 - totalEnergy*totalEnergy);
144  if( nSeeds0 > 6 ) {
145  unsigned nHPD15 = 0;
146  for( const auto& itHPD : theHPDs ) {
147  int hpdN = itHPD.first;
148  const std::vector<unsigned>& hpdHits = itHPD.second;
149  if( ( std::abs(hpdN) < 100 && hpdHits.size() > 14 ) ||
150  ( std::abs(hpdN) > 100 && hpdHits.size() > 14 ) ) ++nHPD15;
151  }
152  if( nHPD15 > 1 ) {
153  unsigned nn = 0;
154  double threshold = 1.0;
155  for( const auto& itEn : theEnergies ) {
156  ++nn;
157  if( nn< 5 ) {
158  mask[itEn.second] = false;
159  } else if ( nn == 5 ) {
160  threshold = itEn.first*5;
161  mask[itEn.second] = false;
162  } else {
163  if( itEn.first < threshold ) mask[itEn.second] = false;
164  }
165  }
166  }
167  }
168  }
169  }
170 
171  // loop on the HPDs
172  std::unordered_map<int,std::vector<unsigned> >::iterator neighbour1;
173  std::unordered_map<int,std::vector<unsigned> >::iterator neighbour2;
174  std::unordered_map<int,std::vector<unsigned> >::iterator neighbour0;
175  std::unordered_map<int,std::vector<unsigned> >::iterator neighbour3;
176  unsigned size1 = 0, size2 = 0;
177  for( const auto& ithpd : _hpds ) {
178  const std::vector<unsigned>& rechits = ithpd.second;
179  theEnergies.clear();
180  totalEnergy = 0;
181  totalEnergy2 = 1e-9;
182  for( const unsigned rhidx : rechits ) {
183  const reco::PFRecHit & rechit = input->at(rhidx);
184  const double e = rechit.energy();
185  totalEnergy += e;
186  totalEnergy2 += e*e;
187  theEnergies.emplace(rechit.energy(),rhidx);
188  }
189  totalEnergy /= rechits.size();
190  totalEnergy2 /= rechits.size();
191  totalEnergy2 = std::sqrt(totalEnergy2 - totalEnergy*totalEnergy);
192 
193  const int thehpd = ithpd.first;
194  switch( std::abs(thehpd) ) {
195  case 1:
196  neighbour1 = ( thehpd > 0 ? _hpds.find(72) : _hpds.find(-72) );
197  break;
198  case 72:
199  neighbour2 = ( thehpd > 0 ? _hpds.find(1) : _hpds.find(-1) );
200  break;
201  case 101:
202  neighbour1 = ( thehpd > 0 ? _hpds.find(136) : _hpds.find(-136) );
203  break;
204  case 136:
205  neighbour2 = ( thehpd > 0 ? _hpds.find(101) : _hpds.find(-101) );
206  break;
207  default:
208  neighbour1 = ( thehpd > 0 ? _hpds.find(thehpd-1) : _hpds.find(thehpd+1) );
209  neighbour2 = ( thehpd > 0 ? _hpds.find(thehpd+1) : _hpds.find(thehpd-1) );
210  break;
211  }
212  if( neighbour1 != _hpds.end() ) {
213  const int nb1 = neighbour1->first;
214  switch( std::abs(nb1) ) {
215  case 1:
216  neighbour0 = ( nb1 > 0 ? _hpds.find(72) : _hpds.find(-72) );
217  break;
218  case 101:
219  neighbour0 = ( nb1 > 0 ? _hpds.find(136) : _hpds.find(-136) );
220  break;
221  default:
222  neighbour0 = ( nb1 > 0 ? _hpds.find(nb1-1) : _hpds.find(nb1+1) );
223  break;
224  }
225  } else {
226  neighbour0 = _hpds.end();
227  }
228 
229  if( neighbour2 != _hpds.end() ) {
230  const int nb2 = neighbour2->first;
231  switch( std::abs(nb2) ) {
232  case 72:
233  neighbour3 = ( nb2 > 0 ? _hpds.find(1) : _hpds.find(-1) );
234  break;
235  case 136:
236  neighbour3 = ( nb2 > 0 ? _hpds.find(101) : _hpds.find(-101) );
237  break;
238  default:
239  neighbour3 = ( nb2 > 0 ? _hpds.find(nb2+1) : _hpds.find(nb2-1) );
240  break;
241  }
242  } else {
243  neighbour3 = _hpds.end();
244  }
245 
246  size1 = neighbour1 != _hpds.end() ? neighbour1->second.size() : 0;
247  size2 = neighbour2 != _hpds.end() ? neighbour2->second.size() : 0;
248  if( size1 > 10 ) {
249  if ( ( abs(neighbour1->first) > 100 && neighbour1->second.size() > 15 ) ||
250  ( abs(neighbour1->first) < 100 && neighbour1->second.size() > 12 ) )
251  size1 = neighbour0 != _hpds.end() ? neighbour0->second.size() : 0;
252  }
253  if( size2 > 10 ) {
254  if ( ( abs(neighbour2->first) > 100 && neighbour2->second.size() > 15 ) ||
255  ( abs(neighbour2->first) < 100 && neighbour2->second.size() > 12 ) )
256  size2 = neighbour3 != _hpds.end() ? neighbour3->second.size() : 0;
257  }
258  if( ( std::abs(ithpd.first) > 100 && ithpd.second.size() > 15 ) ||
259  ( std::abs(ithpd.first) < 100 && ithpd.second.size() > 12 ) ) {
260  if( (double)(size1+size2)/(float)ithpd.second.size() < 1.0 ) {
261  unsigned nn = 0;
262  double threshold = 1.0;
263  for( const auto& itEn : theEnergies ) {
264  if( nn < 5 ) {
265  mask[itEn.second] = false;
266  } else if ( nn == 5 ) {
267  threshold = itEn.first*2.5;
268  mask[itEn.second] = false;
269  } else {
270  if( itEn.first < threshold ) mask[itEn.second] = false;
271  }
272  }
273  }
274  }
275  }
276 }
i
int i
Definition: DBlmapReader.cc:9
reco::PFRecHit::neighbours4
const PFRecHitRefVector & neighbours4() const
Definition: PFRecHit.h:79
reco::PFRecHit::detId
unsigned detId() const
rechit detId
Definition: PFRecHit.h:101
reco::PFRecHit::layer
PFLayer::Layer layer() const
rechit layer
Definition: PFRecHit.h:104
reco::PFRecHit
Particle flow rechit (rechit + geometry and topology information). See clustering algorithm in PFClus...
Definition: PFRecHit.h:35
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:48
HcalDetId::ieta
int ieta() const
get the cell ieta
Definition: HcalDetId.h:36
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
dtDQMClient_cfg.threshold
tuple threshold
Definition: dtDQMClient_cfg.py:16
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.h:38
reco::PFRecHit::position
const math::XYZPoint & position() const
rechit cell centre x, y, z
Definition: PFRecHit.h:123
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::energy
double energy() const
rechit energy
Definition: PFRecHit.h:107
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