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

#include <PreshowerClusterAlgo.h>

Public Types

typedef std::set< DetId > HitsID
 
typedef math::XYZPoint Point
 
typedef std::map< DetId,
EcalRecHit > 		RecHitsMap
 

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_
 
RecHitsMap * rechits_map
 
std::vector< ESDetId > road_2d
 
HitsID * used_s
 

Detailed Description

Definition at line 23 of file PreshowerClusterAlgo.h.

Member Typedef Documentation

typedef std::set<DetId> PreshowerClusterAlgo::HitsID

Definition at line 32 of file PreshowerClusterAlgo.h.

typedef math::XYZPoint PreshowerClusterAlgo::Point

Definition at line 29 of file PreshowerClusterAlgo.h.

typedef std::map<DetId, EcalRecHit> PreshowerClusterAlgo::RecHitsMap

Definition at line 31 of file PreshowerClusterAlgo.h.

Constructor & Destructor Documentation

PreshowerClusterAlgo::PreshowerClusterAlgo ( )
inline

Definition at line 34 of file PreshowerClusterAlgo.h.

34  :
35  preshStripEnergyCut_(0.), preshClusterEnergyCut_(0.), preshSeededNstr_(15)
36  {}
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.

42 {};

Member Function Documentation

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

Definition at line 268 of file PreshowerClusterAlgo.cc.

References LogTrace, GetRecoTauVFromDQM_MC_cff::next, preshSeededNstr_, and road_2d.

Referenced by makeOneCluster().

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

Definition at line 246 of file PreshowerClusterAlgo.cc.

References LogTrace, preshStripEnergyCut_, rechits_map, and used_s.

Referenced by makeOneCluster().

247 {
248 
249  if ( used_s->find(candidate_it->first) != used_s->end())
250  LogTrace("PreShowerClusterAlgo") << " This strip is in use";
251  if (candidate_it == rechits_map->end() )
252  LogTrace("PreShowerClusterAlgo") << " No such a strip in rechits_map";
253  if (candidate_it->second.energy() <= preshStripEnergyCut_)
254  LogTrace("PreShowerClusterAlgo") << "Strip energy" << candidate_it->second.energy() <<"is below threshold";
255 
256  // crystal should not be included...
257  if ( (used_s->find(candidate_it->first) != used_s->end()) || // ...if it already belongs to a cluster
258  (candidate_it == rechits_map->end() ) || // ...if it corresponds to a hit
259  (candidate_it->second.energy() <= preshStripEnergyCut_ ) ) // ...if it has a negative or zero energy
260  {
261  return false;
262  }
263  return true;
264 }
LogTrace
#define LogTrace(id)
Definition: MessageLogger.h:502
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(), CommonMethods::cp(), edm::SortedCollection< T, SORT >::end(), reco::CaloCluster::energy(), reco::CaloCluster::eta(), findRoad(), newFWLiteAna::found, CaloSubdetectorGeometry::getGeometry(), CaloCellGeometry::getPosition(), goodStrip(), LogTrace, GetRecoTauVFromDQM_MC_cff::next, reco::PreshowerCluster::nhits(), reco::CaloCluster::phi(), ESDetId::plane(), reco::CaloCluster::position(), position, preshClusterEnergyCut_, edm::SortedCollection< T, SORT >::push_back(), rechits_map, road_2d, edm::SortedCollection< T, SORT >::size(), ESDetId::strip(), used_s, 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().

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  const CaloCellGeometry *thisCell = geometry_p->getGeometry(cp->first);
201  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 
214  EcalRecHitCollection::iterator it;
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  used_strips = used_s;
236 
237  // return the cluster if its energy is greater a threshold
238  if( cluster.energy() > preshClusterEnergyCut_ )
239  return cluster;
240  else
241  return nullcluster;
242 
243 }
EcalPreshowerNavigator
CaloNavigator< ESDetId > EcalPreshowerNavigator
Definition: EcalPreshowerNavigator.h:7
reco::CaloCluster::position
const math::XYZPoint & position() const
cluster centroid position
Definition: CaloCluster.h:124
ESDetId::strip
int strip() const
Definition: ESDetId.h:52
PreshowerClusterAlgo::findRoad
void findRoad(ESDetId strip, EcalPreshowerNavigator theESNav, int plane)
Definition: PreshowerClusterAlgo.cc:268
PV3DBase::y
T y() const
Definition: PV3DBase.h:63
edm::SortedCollection::push_back
void push_back(T const &t)
Definition: SortedCollection.h:190
reco::CaloCluster::z
double z() const
z coordinate of cluster centroid
Definition: CaloCluster.h:158
position
static int position[TOTALCHAMBERS][3]
Definition: ReadPGInfo.cc:509
reco::CaloCluster::eta
double eta() const
pseudorapidity of cluster centroid
Definition: CaloCluster.h:161
CaloSubdetectorGeometry::getGeometry
virtual const CaloCellGeometry * getGeometry(const DetId &id) const
Get the cell geometry of a given detector id. Should return false if not found.
Definition: CaloSubdetectorGeometry.cc:44
RecHitsMap
std::map< DetId, EcalRecHit > RecHitsMap
Definition: HLTAlCaMonPi0.h:24
PV3DBase::z
T z() const
Definition: PV3DBase.h:64
Point
math::XYZPoint Point
Definition: TrackerDpgAnalysis.cc:107
PreshowerClusterAlgo::goodStrip
bool goodStrip(RecHitsMap::iterator candidate_it)
Definition: PreshowerClusterAlgo.cc:246
reco::CaloCluster::x
double x() const
x coordinate of cluster centroid
Definition: CaloCluster.h:152
reco::CaloCluster::energy
double energy() const
cluster energy
Definition: CaloCluster.h:120
reco::PreshowerCluster::nhits
int nhits() const
Number of RecHits the cluster.
Definition: PreshowerCluster.h:35
LogTrace
#define LogTrace(id)
Definition: MessageLogger.h:502
PreshowerClusterAlgo::road_2d
std::vector< ESDetId > road_2d
Definition: PreshowerClusterAlgo.h:61
edm::SortedCollection< EcalRecHit >::iterator
std::vector< EcalRecHit >::iterator iterator
Definition: SortedCollection.h:81
edm::SortedCollection::end
const_iterator end() const
Definition: SortedCollection.h:301
edm::SortedCollection::size
size_type size() const
Definition: SortedCollection.h:227
ESDetId::plane
int plane() const
Definition: ESDetId.h:46
GetRecoTauVFromDQM_MC_cff.next
tuple next
Definition: GetRecoTauVFromDQM_MC_cff.py:30
reco::CaloCluster::y
double y() const
y coordinate of cluster centroid
Definition: CaloCluster.h:155
reco::CaloCluster::phi
double phi() const
azimuthal angle of cluster centroid
Definition: CaloCluster.h:164
CaloCellGeometry::getPosition
const GlobalPoint & getPosition() const
Returns the position of reference for this cell.
Definition: CaloCellGeometry.h:76
PV3DBase::x
T x() const
Definition: PV3DBase.h:62
edm::SortedCollection::begin
const_iterator begin() const
Definition: SortedCollection.h:294

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(), and makeOneCluster().

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().

Generated for CMSSW Reference Manual by doxygen 1.8.5