CMS 3D CMS Logo

List of all members | Public Types | Public Member Functions | Private Attributes
PreshowerClusterAlgo Class Reference

#include <PreshowerClusterAlgo.h>

Public Types

typedef std::set< DetIdHitsID
 
typedef math::XYZPoint Point
 
typedef std::map< DetId, EcalRecHitRecHitsMap
 

Public Member Functions

void findRoad (ESDetId strip, EcalPreshowerNavigator theESNav, int plane)
 
bool goodStrip (RecHitsMap::iterator candidate_it)
 
reco::PreshowerCluster makeOneCluster (ESDetId strip, HitsID *used_strips, RecHitsMap *rechits_map, const CaloSubdetectorGeometry *&geometry_p, CaloSubdetectorTopology *&topology_p)
 
 PreshowerClusterAlgo ()
 
 PreshowerClusterAlgo (double stripEnergyCut, double clusterEnergyCut, int nStripCut)
 
 ~PreshowerClusterAlgo ()
 

Private Attributes

double preshClusterEnergyCut_
 
int preshSeededNstr_
 
double preshStripEnergyCut_
 
RecHitsMaprechits_map
 
std::vector< ESDetIdroad_2d
 
HitsIDused_s
 

Detailed Description

Definition at line 23 of file PreshowerClusterAlgo.h.

Member Typedef Documentation

Definition at line 32 of file PreshowerClusterAlgo.h.

Definition at line 29 of file PreshowerClusterAlgo.h.

Definition at line 31 of file PreshowerClusterAlgo.h.

Constructor & Destructor Documentation

PreshowerClusterAlgo::PreshowerClusterAlgo ( )
inline
PreshowerClusterAlgo::PreshowerClusterAlgo ( double  stripEnergyCut,
double  clusterEnergyCut,
int  nStripCut 
)
inline

Definition at line 38 of file PreshowerClusterAlgo.h.

38  :
39  preshStripEnergyCut_(stripEnergyCut), preshClusterEnergyCut_(clusterEnergyCut), preshSeededNstr_(nStripCut)
40  {}
PreshowerClusterAlgo::~PreshowerClusterAlgo ( )
inline

Definition at line 42 of file PreshowerClusterAlgo.h.

References findRoad(), goodStrip(), makeOneCluster(), rechits_map, and digitizers_cfi::strip.

42 {};

Member Function Documentation

void PreshowerClusterAlgo::findRoad ( ESDetId  strip,
EcalPreshowerNavigator  theESNav,
int  plane 
)

Definition at line 267 of file PreshowerClusterAlgo.cc.

References CaloNavigator< T, TOPO >::east(), CaloNavigator< T, TOPO >::home(), LogTrace, GetRecoTauVFromDQM_MC_cff::next, CaloNavigator< T, TOPO >::north(), preshSeededNstr_, road_2d, CaloNavigator< T, TOPO >::setHome(), CaloNavigator< T, TOPO >::south(), digitizers_cfi::strip, and CaloNavigator< T, TOPO >::west().

Referenced by makeOneCluster(), and ~PreshowerClusterAlgo().

267  {
268 
269  if ( strip == ESDetId(0) ) return;
270 
271  ESDetId next;
272  theESNav.setHome(strip);
273 // First, add a central strip to the road
274  road_2d.push_back(strip);
275  LogTrace("PreShowerClusterAlgo") << "findRoad starts from strip" << strip;
276 
277  if (plane == 1) {
278  // east road
279  int n_east= 0;
280  LogTrace("PreShowerClusterAlgo") << " Go to the East ";
281 
282  while ( ((next=theESNav.east()) != ESDetId(0) && next != strip) ) {
283  LogTrace("PreShowerClusterAlgo") << "East:" << n_east << "current strip is"<< next;
284 
285  road_2d.push_back(next);
286  ++n_east;
287  if (n_east == preshSeededNstr_) break;
288  }
289  // west road
290  int n_west= 0;
291  LogTrace("PreShowerClusterAlgo") << "Go to the West" ;
292 
293  theESNav.home();
294  while ( ((next=theESNav.west()) != ESDetId(0) && next != strip )) {
295  LogTrace("PreShowerClusterAlgo") << "West: " << n_west << "current strip is" << next ;
296 
297  road_2d.push_back(next);
298  ++n_west;
299  if (n_west == preshSeededNstr_) break;
300  }
301  LogTrace("PreShowerClusterAlgo") << "Total number of strips found in the road at 1-st plane is" << n_east+n_west ;
302 
303  }
304  else if (plane == 2) {
305  // north road
306  int n_north= 0;
307  LogTrace("PreShowerClusterAlgo") << "Go to the North";
308 
309  while ( ((next=theESNav.north()) != ESDetId(0) && next != strip) ) {
310  LogTrace("PreShowerClusterAlgo") << "North:" << n_north << "current strip is" << next ;
311 
312  road_2d.push_back(next);
313  ++n_north;
314  if (n_north == preshSeededNstr_) break;
315  }
316  // south road
317  int n_south= 0;
318  LogTrace("PreShowerClusterAlgo") << "Go to the South";
319 
320  theESNav.home();
321  while ( ((next=theESNav.south()) != ESDetId(0) && next != strip) ) {
322 
323  LogTrace("PreShowerClusterAlgo") << "South:" << n_south << "current strip is" << next ;
324 
325  road_2d.push_back(next);
326  ++n_south;
327  if (n_south == preshSeededNstr_) break;
328  }
329 
330  LogTrace("PreShowerClusterAlgo") << "Total number of strips found in the road at 2-nd plane is" << n_south+n_north;
331 
332  }
333  else {
334  LogTrace("PreShowerClusterAlgo") << " Wrong plane number, null cluster will be returned!";
335 
336  } // end of if
337 
338  theESNav.home();
339 }
void setHome(const T &startingPoint)
set the starting position
T west() const
move the navigator west
Definition: CaloNavigator.h:59
T south() const
move the navigator south
Definition: CaloNavigator.h:45
#define LogTrace(id)
std::vector< ESDetId > road_2d
T east() const
move the navigator east
Definition: CaloNavigator.h:52
void home() const
move the navigator back to the starting point
T north() const
move the navigator north
Definition: CaloNavigator.h:38
bool PreshowerClusterAlgo::goodStrip ( RecHitsMap::iterator  candidate_it)

Definition at line 245 of file PreshowerClusterAlgo.cc.

References LogTrace, preshStripEnergyCut_, rechits_map, and used_s.

Referenced by makeOneCluster(), and ~PreshowerClusterAlgo().

246 {
247 
248  if ( used_s->find(candidate_it->first) != used_s->end())
249  LogTrace("PreShowerClusterAlgo") << " This strip is in use";
250  if (candidate_it == rechits_map->end() )
251  LogTrace("PreShowerClusterAlgo") << " No such a strip in rechits_map";
252  if (candidate_it->second.energy() <= preshStripEnergyCut_)
253  LogTrace("PreShowerClusterAlgo") << "Strip energy" << candidate_it->second.energy() <<"is below threshold";
254 
255  // crystal should not be included...
256  if ( (used_s->find(candidate_it->first) != used_s->end()) || // ...if it already belongs to a cluster
257  (candidate_it == rechits_map->end() ) || // ...if it corresponds to a hit
258  (candidate_it->second.energy() <= preshStripEnergyCut_ ) ) // ...if it has a negative or zero energy
259  {
260  return false;
261  }
262  return true;
263 }
#define LogTrace(id)
reco::PreshowerCluster PreshowerClusterAlgo::makeOneCluster ( ESDetId  strip,
HitsID used_strips,
RecHitsMap rechits_map,
const CaloSubdetectorGeometry *&  geometry_p,
CaloSubdetectorTopology *&  topology_p 
)

Definition at line 13 of file PreshowerClusterAlgo.cc.

References edm::SortedCollection< T, SORT >::begin(), SimDataFormats::CaloAnalysis::cp, CaloNavigator< T, TOPO >::east(), edm::SortedCollection< T, SORT >::end(), CaloParticle::energy(), reco::CaloCluster::energy(), reco::CaloCluster::eta(), findRoad(), runEdmFileComparison::found, CaloSubdetectorGeometry::getGeometry(), goodStrip(), CaloNavigator< T, TOPO >::home(), LogTrace, particleFlowRecHitECAL_cfi::navigator, GetRecoTauVFromDQM_MC_cff::next, reco::PreshowerCluster::nhits(), CaloNavigator< T, TOPO >::north(), reco::CaloCluster::phi(), ESDetId::plane(), reco::CaloCluster::position(), position, preshClusterEnergyCut_, edm::SortedCollection< T, SORT >::push_back(), rechits_map, road_2d, CaloNavigator< T, TOPO >::setHome(), edm::SortedCollection< T, SORT >::size(), CaloNavigator< T, TOPO >::south(), ESDetId::strip(), used_s, CaloNavigator< T, TOPO >::west(), PV3DBase< T, PVType, FrameType >::x(), reco::CaloCluster::x(), PV3DBase< T, PVType, FrameType >::y(), reco::CaloCluster::y(), PV3DBase< T, PVType, FrameType >::z(), and reco::CaloCluster::z().

Referenced by ~PreshowerClusterAlgo().

18 {
19  road_2d.clear();
20 
21  rechits_map = the_rechitsMap_p;
22 
23  used_s = used_strips;
24 
25  int plane = strip.plane();
26 
27  LogTrace("PreShowerClusterAlgo") << "Preshower Seeded Algorithm - looking for clusters";
28  LogTrace("PreShowerClusterAlgo")<< "Preshower is intersected at strip" << strip.strip() << ",at plane" << plane ;
29 
30 
31  // create null-cluster
32  std::vector< std::pair<DetId,float> > dummy;
33  Point posi(0,0,0);
34  LogTrace("PreShowerClusterAlgo") << " Creating a null-cluster" ;
35 
36  reco::PreshowerCluster nullcluster=reco::PreshowerCluster(0.,posi,dummy,plane);
37 
38  if ( strip == ESDetId(0) ) return nullcluster; //works in case of no intersected strip found (e.g. in the Barrel)
39 
40  // Collection of cluster strips
41  EcalRecHitCollection clusterRecHits;
42  // Map of strips for position calculation
43  RecHitsMap recHits_pos;
44 
45  //Make a navigator, and set it to the strip cell.
46  EcalPreshowerNavigator navigator(strip, topology_p);
47  navigator.setHome(strip);
48  //search for neighbours in the central road
49  findRoad(strip,navigator,plane);
50  LogTrace("PreShowerClusterAlgo") << "Total number of strips in the central road:" << road_2d.size() ;
51 
52  if ( plane == 1 ) {
53  ESDetId strip_north = navigator.north();
54  findRoad(strip_north,navigator,plane);
55  navigator.home();
56  ESDetId strip_south = navigator.south();
57  findRoad(strip_south,navigator,plane);
58  navigator.home();
59  }
60  if ( plane == 2 ) {
61  ESDetId strip_east = navigator.east();
62  findRoad(strip_east,navigator,plane);
63  navigator.home();
64  ESDetId strip_west = navigator.west();
65  findRoad(strip_west,navigator,plane);
66  navigator.home();
67  }
68  LogTrace("PreShowerClusterAlgo") << "Total number of strips in all three roads:" << road_2d.size() ;
69 
70  // Start clustering from strip with max Energy in the road
71  float E_max = 0.;
72  bool found = false;
73  RecHitsMap::iterator max_it;
74  // Loop over strips:
75  std::vector<ESDetId>::iterator itID;
76  for (itID = road_2d.begin(); itID != road_2d.end(); itID++) {
77  LogTrace("PreShowerClusterAlgo") << "ID ="<<*itID ;
78 
79  RecHitsMap::iterator strip_it = rechits_map->find(*itID);
80  //if ( strip_it->second.energy() < 0 ) std::cout << " ##### E = " << strip_it->second.energy() << std::endl;
81  if( strip_it==rechits_map->end() || !goodStrip(strip_it)) continue;
82  LogTrace("PreShowerClusterAlgo") << " strip is " << ESDetId(strip_it->first) <<"E ="<< strip_it->second.energy();
83 
84  float E = strip_it->second.energy();
85  if ( E > E_max) {
86  E_max = E;
87  found = true;
88  max_it = strip_it;
89  }
90  }
91  // no hottest strip found ==> null cluster will be returned
92  if ( !found ) return nullcluster;
93 
94  // First, save the hottest strip
95  clusterRecHits.push_back(max_it->second);
96  recHits_pos.insert(std::make_pair(max_it->first, max_it->second));
97  used_s->insert(max_it->first);
98  LogTrace("PreShowerClusterAlgo") << "Central hottest strip" << ESDetId(max_it->first) << "is saved with energy E =" << E_max ;
99 
100 
101  // Find positions of adjacent strips:
102  ESDetId next, strip_1, strip_2;
103  navigator.setHome(max_it->first);
104  ESDetId startES = max_it->first;
105 
106  if (plane == 1) {
107  // Save two neighbouring strips to the east
108  int nadjacents_east = 0;
109  while ( (next=navigator.east()) != ESDetId(0) && next != startES && nadjacents_east < 2 ) {
110  ++nadjacents_east;
111  LogTrace("PreShowerClusterAlgo") << "Adjacent east #" << nadjacents_east <<":"<< next ;
112 
113  RecHitsMap::iterator strip_it = rechits_map->find(next);
114 
115  if(strip_it==rechits_map->end() || !goodStrip(strip_it)) continue;
116  // Save strip for clustering if it exists, not already in use, and satisfies an energy threshold
117  clusterRecHits.push_back(strip_it->second);
118  // save strip for position calculation
119  if ( nadjacents_east==1 ) strip_1 = next;
120  used_s->insert(strip_it->first);
121  LogTrace("PreShowerClusterAlgo") << "East adjacent strip #" << nadjacents_east << "is saved with energy E =" << strip_it->second.energy() ;
122 
123  }
124  // Save two neighbouring strips to the west
125  navigator.home();
126  int nadjacents_west = 0;
127  while ( (next=navigator.west()) != ESDetId(0) && next != startES && nadjacents_west < 2 ) {
128  ++nadjacents_west;
129  LogTrace("PreShowerClusterAlgo") << "Adjacent west #" << nadjacents_west <<":"<< next ;
130 
131  RecHitsMap::iterator strip_it = rechits_map->find(next);
132  if(strip_it==rechits_map->end() || !goodStrip(strip_it)) continue;
133  clusterRecHits.push_back(strip_it->second);
134  if ( nadjacents_west==1 ) strip_2 = next;
135  used_s->insert(strip_it->first);
136  LogTrace("PreShowerClusterAlgo") << "West adjacent strip #" << nadjacents_west << "is saved with energy E =" << strip_it->second.energy();
137 
138  }
139  }
140  else if (plane == 2) {
141 
142  // Save two neighbouring strips to the north
143  int nadjacents_north = 0;
144  while ( (next=navigator.north()) != ESDetId(0) && next != startES && nadjacents_north < 2 ) {
145  ++nadjacents_north;
146  LogTrace("PreShowerClusterAlgo") << "Adjacent north #" << nadjacents_north <<":"<< next ;
147 
148  RecHitsMap::iterator strip_it = rechits_map->find(next);
149  if(strip_it==rechits_map->end() || !goodStrip(strip_it)) continue;
150  clusterRecHits.push_back(strip_it->second);
151  if ( nadjacents_north==1 ) strip_1 = next;
152  used_s->insert(strip_it->first);
153  LogTrace("PreShowerClusterAlgo") << "North adjacent strip #" << nadjacents_north << "is saved with energy E =" << strip_it->second.energy() ;
154 
155  }
156  // Save two neighbouring strips to the south
157  navigator.home();
158  int nadjacents_south = 0;
159  while ( (next=navigator.south()) != ESDetId(0) && next != startES && nadjacents_south < 2 ) {
160  ++nadjacents_south;
161  LogTrace("PreShowerClusterAlgo") << "Adjacent south #" << nadjacents_south <<":"<< next ;
162 
163  RecHitsMap::iterator strip_it = rechits_map->find(next);
164  if(strip_it==rechits_map->end() || !goodStrip(strip_it)) continue;
165  clusterRecHits.push_back(strip_it->second);
166  if ( nadjacents_south==1 ) strip_2 = next;
167  used_s->insert(strip_it->first);
168  LogTrace("PreShowerClusterAlgo") << "South adjacent strip #" << nadjacents_south << "is saved with energy E =" << strip_it->second.energy() ;
169 
170  }
171  }
172  else {
173  LogTrace("PreShowerClusterAlgo") << " Wrong plane number" << plane <<", null cluster will be returned! " << std::endl;
174  return nullcluster;
175  } // end of if
176 
177  LogTrace("PreShowerClusterAlgo") << "Total size of clusterRecHits is" << clusterRecHits.size();
178  LogTrace("PreShowerClusterAlgo") << "Two adjacent strips for position calculation are:" << strip_1 <<"and" << strip_2;
179 
180 
181  // strips for position calculation
182  RecHitsMap::iterator strip_it1, strip_it2;
183  if ( strip_1 != ESDetId(0)) {
184  strip_it1 = rechits_map->find(strip_1);
185  recHits_pos.insert(std::make_pair(strip_it1->first, strip_it1->second));
186  }
187  if ( strip_2 != ESDetId(0) ) {
188  strip_it2 = rechits_map->find(strip_2);
189  recHits_pos.insert(std::make_pair(strip_it2->first, strip_it2->second));
190  }
191 
192  RecHitsMap::iterator cp;
193  double energy_pos = 0;
194  double x_pos = 0;
195  double y_pos = 0;
196  double z_pos = 0;
197  for (cp = recHits_pos.begin(); cp!=recHits_pos.end(); cp++ ) {
198  double E = cp->second.energy();
199  energy_pos += E;
200  auto thisCell = geometry_p->getGeometry(cp->first);
201  const GlobalPoint& position = thisCell->getPosition();
202  x_pos += E * position.x();
203  y_pos += E * position.y();
204  z_pos += E * position.z();
205  }
206  if(energy_pos>0.) {
207  x_pos /= energy_pos;
208  y_pos /= energy_pos;
209  z_pos /= energy_pos;
210  }
211  Point pos(x_pos,y_pos,z_pos);
212  LogTrace("PreShowerClusterAlgo") << "ES Cluster position =" << "("<< x_pos <<","<< y_pos <<","<< z_pos <<")";
213 
215  double Eclust = 0;
216 
217 
218  std::vector<std::pair<DetId,float > > usedHits;
219  for (it=clusterRecHits.begin(); it != clusterRecHits.end(); it++) {
220  Eclust += it->energy();
221  usedHits.push_back(std::pair<DetId,float > (it->id(),1.));
222 
223  }
224 
225 
226  // ES cluster is created from vector clusterRecHits
227  reco::PreshowerCluster cluster=reco::PreshowerCluster(Eclust,pos,usedHits,plane);
228  LogTrace("PreShowerClusterAlgo") << " ES Cluster is created with:" ;
229  LogTrace("PreShowerClusterAlgo") << " energy =" << cluster.energy();
230  LogTrace("PreShowerClusterAlgo") << " (eta,phi) =" << "("<<cluster.eta()<<","<<cluster.phi()<<")";
231  LogTrace("PreShowerClusterAlgo") << " nhits =" << cluster.nhits();
232  LogTrace("PreShowerClusterAlgo") << " radius =" << cluster.position().r();
233  LogTrace("PreShowerClusterAlgo") << " (x,y,z) =" << "(" << cluster.x() <<", "<< cluster.y() <<","<< cluster.z()<<")" ;
234 
235 
236  // return the cluster if its energy is greater a threshold
237  if( cluster.energy() > preshClusterEnergyCut_ )
238  return cluster;
239  else
240  return nullcluster;
241 
242 }
const math::XYZPoint & position() const
cluster centroid position
Definition: CaloCluster.h:131
int strip() const
Definition: ESDetId.h:53
void findRoad(ESDetId strip, EcalPreshowerNavigator theESNav, int plane)
T y() const
Definition: PV3DBase.h:63
void push_back(T const &t)
double z() const
z coordinate of cluster centroid
Definition: CaloCluster.h:165
double eta() const
pseudorapidity of cluster centroid
Definition: CaloCluster.h:168
T z() const
Definition: PV3DBase.h:64
math::XYZPoint Point
float energy() const
Energy. Note this is taken from the first SimTrack only.
Definition: CaloParticle.h:106
bool goodStrip(RecHitsMap::iterator candidate_it)
double x() const
x coordinate of cluster centroid
Definition: CaloCluster.h:159
double energy() const
cluster energy
Definition: CaloCluster.h:126
int nhits() const
Number of RecHits the cluster.
#define LogTrace(id)
std::vector< ESDetId > road_2d
std::vector< EcalRecHit >::iterator iterator
const_iterator end() const
virtual std::shared_ptr< const CaloCellGeometry > getGeometry(const DetId &id) const
Get the cell geometry of a given detector id. Should return false if not found.
size_type size() const
static int position[264][3]
Definition: ReadPGInfo.cc:509
int plane() const
Definition: ESDetId.h:47
double y() const
y coordinate of cluster centroid
Definition: CaloCluster.h:162
double phi() const
azimuthal angle of cluster centroid
Definition: CaloCluster.h:171
T x() const
Definition: PV3DBase.h:62
const_iterator begin() const
std::map< DetId, EcalRecHit > RecHitsMap
Definition: DQMSourcePi0.h:27

Member Data Documentation

double PreshowerClusterAlgo::preshClusterEnergyCut_
private

Definition at line 57 of file PreshowerClusterAlgo.h.

Referenced by makeOneCluster().

int PreshowerClusterAlgo::preshSeededNstr_
private

Definition at line 58 of file PreshowerClusterAlgo.h.

Referenced by findRoad().

double PreshowerClusterAlgo::preshStripEnergyCut_
private

Definition at line 56 of file PreshowerClusterAlgo.h.

Referenced by goodStrip().

RecHitsMap* PreshowerClusterAlgo::rechits_map
private

Definition at line 64 of file PreshowerClusterAlgo.h.

Referenced by goodStrip(), makeOneCluster(), and ~PreshowerClusterAlgo().

std::vector<ESDetId> PreshowerClusterAlgo::road_2d
private

Definition at line 61 of file PreshowerClusterAlgo.h.

Referenced by findRoad(), and makeOneCluster().

HitsID* PreshowerClusterAlgo::used_s
private

Definition at line 67 of file PreshowerClusterAlgo.h.

Referenced by goodStrip(), and makeOneCluster().