CMS 3D CMS Logo

Multi5x5ClusterAlgo.cc
Go to the documentation of this file.
1 
3 
4 // Geometry
11 
14 
15 //temporary sorter, I'm sure this must exist already in CMSSW
16 template <class T1,class T2,typename Comp=std::less<T1> > struct PairSortByFirst {
17  Comp comp;
18  bool operator()(const std::pair<T1,T2>& lhs,const std::pair<T1,T2>&rhs){return comp(lhs.first,rhs.first);}
19  bool operator()(const T1& lhs,const std::pair<T1,T2>&rhs){return comp(lhs,rhs.first);}
20  bool operator()(const std::pair<T1,T2>& lhs,const T1 &rhs){return comp(lhs.first,rhs);}
21  bool operator()(const T1& lhs,const T1 &rhs){return comp(lhs,rhs);}
22 };
23 
25  bool operator()(const EBDetId& lhs,const EBDetId& rhs){
26  if(lhs.ieta()<rhs.ieta()) return true;
27  else if(lhs.ieta()>rhs.ieta()) return false;
28  else return lhs.iphi()<rhs.iphi();
29  }
30 };
31 
33  bool operator()(const EEDetId& lhs,const EEDetId& rhs){
34  if(lhs.zside()<rhs.zside()) return true;
35  else if(lhs.zside()>rhs.zside()) return false;
36  else { //z is equal, onto ix
37  if(lhs.ix()<rhs.ix()) return true;
38  else if(lhs.ix()>rhs.ix()) return false;
39  else return lhs.iy()<rhs.iy();
40  }
41  }
42 };
43 
44 // Return a vector of clusters from a collection of EcalRecHits:
45 //
46 std::vector<reco::BasicCluster> Multi5x5ClusterAlgo::makeClusters(
48  const CaloSubdetectorGeometry *geometry_p,
49  const CaloSubdetectorTopology *topology_p,
50  const CaloSubdetectorGeometry *geometryES_p,
52  bool regional,
53  const std::vector<EcalEtaPhiRegion>& regions
54  )
55 {
56  seeds.clear();
57  used_s.clear();
58  canSeed_s.clear();
59  clusters_v.clear();
60 
61  recHits_ = hits;
62 
63  double threshold = 0;
64  std::string ecalPart_string;
65  detector_ = reco::CaloID::DET_NONE;
66  if (detector == reco::CaloID::DET_ECAL_ENDCAP)
67  {
69  threshold = ecalEndcapSeedThreshold;
70  ecalPart_string = "EndCap";
71  }
72  if (detector == reco::CaloID::DET_ECAL_BARREL)
73  {
75  threshold = ecalBarrelSeedThreshold;
76  ecalPart_string = "Barrel";
77  }
78 
79  LogTrace("EcalClusters") << "-------------------------------------------------------------";
80  LogTrace("EcalClusters") << "Island algorithm invoked for ECAL" << ecalPart_string ;
81  LogTrace("EcalClusters") << "Looking for seeds, energy threshold used = " << threshold << " GeV";
82 
83 
84  int nregions=0;
85  if(regional) nregions=regions.size();
86 
87  if(!regional || nregions) {
88 
90  for(it = hits->begin(); it != hits->end(); it++)
91  {
92  double energy = it->energy();
93  if (energy < threshold) continue; // need to check to see if this line is useful!
94 
95  auto thisCell = geometry_p->getGeometry(it->id());
96  // Require that RecHit is within clustering region in case
97  // of regional reconstruction
98  bool withinRegion = false;
99  if (regional) {
100  std::vector<EcalEtaPhiRegion>::const_iterator region;
101  for (region=regions.begin(); region!=regions.end(); region++) {
102  if (region->inRegion(thisCell->etaPos(),thisCell->phiPos())) {
103  withinRegion = true;
104  break;
105  }
106  }
107  }
108 
109  if (!regional || withinRegion) {
110  float ET = it->energy() * thisCell->getPosition().basicVector().unit().perp();
111  if (ET > threshold) seeds.push_back(*it);
112  }
113  }
114 
115  }
116 
117  sort(seeds.begin(), seeds.end(), EcalRecHitLess());
118 
119  LogTrace("EcalClusters") << "Total number of seeds found in event = " << seeds.size();
120 
121 
122  mainSearch(hits, geometry_p, topology_p, geometryES_p);
123  sort(clusters_v.rbegin(), clusters_v.rend(), ClusterEtLess());
124 
125  LogTrace("EcalClusters") << "---------- end of main search. clusters have been sorted ----";
126 
127 
128  return clusters_v;
129 
130 }
131 
132 // Search for clusters
133 //
135  const CaloSubdetectorGeometry *geometry_p,
136  const CaloSubdetectorTopology *topology_p,
137  const CaloSubdetectorGeometry *geometryES_p
138  )
139 {
140 
141  LogTrace("EcalClusters") << "Building clusters............";
142 
143  // Loop over seeds:
144  std::vector<EcalRecHit>::iterator it;
145  for (it = seeds.begin(); it != seeds.end(); it++)
146  {
147 
148  // check if this crystal is able to seed
149  // (event though it is already used)
150  bool usedButCanSeed = false;
151  if (canSeed_s.find(it->id()) != canSeed_s.end()) usedButCanSeed = true;
152 
153  // avoid seeding for anomalous channels
154  if(! it->checkFlag(EcalRecHit::kGood)){ // if rechit is good, no need for further checks
155  if (it->checkFlags( v_chstatus_ )) continue;
156  }
157 
158  // make sure the current seed does not belong to a cluster already.
159  if ((used_s.find(it->id()) != used_s.end()) && (usedButCanSeed == false))
160  {
161  if (it == seeds.begin())
162  {
163  LogTrace("EcalClusters") << "##############################################################" ;
164  LogTrace("EcalClusters") << "DEBUG ALERT: Highest energy seed already belongs to a cluster!";
165  LogTrace("EcalClusters") << "##############################################################";
166 
167  }
168 
169  // seed crystal is used or is used and cannot seed a cluster
170  // so continue to the next seed crystal...
171  continue;
172  }
173 
174  // clear the vector of hits in current cluster
175  current_v.clear();
176 
177  // Create a navigator at the seed and get seed
178  // energy
179  CaloNavigator<DetId> navigator(it->id(), topology_p);
180  DetId seedId = navigator.pos();
181  EcalRecHitCollection::const_iterator seedIt = hits->find(seedId);
182  navigator.setHome(seedId);
183 
184  // Is the seed a local maximum?
185  bool localMaxima = checkMaxima(navigator, hits);
186 
187  if (localMaxima)
188  {
189  // build the 5x5 taking care over which crystals
190  // can seed new clusters and which can't
191  prepareCluster(navigator, hits, geometry_p);
192  }
193 
194  // If some crystals in the current vector then
195  // make them into a cluster
196  if (!current_v.empty())
197  {
198  makeCluster(hits, geometry_p, geometryES_p, seedIt, usedButCanSeed);
199  }
200 
201  } // End loop on seed crystals
202 
203 
204  if(reassignSeedCrysToClusterItSeeds_){
205  std::sort(whichClusCrysBelongsTo_.begin(),whichClusCrysBelongsTo_.end(),PairSortByFirst<DetId,int>());
206 
207 
208  for(size_t clusNr=0;clusNr<protoClusters_.size();clusNr++){
209  if(!protoClusters_[clusNr].containsSeed()){
210  const EcalRecHit& seedHit =protoClusters_[clusNr].seed();
211  typedef std::vector<std::pair<DetId,int> >::iterator It;
212  std::pair<It,It> result = std::equal_range(whichClusCrysBelongsTo_.begin(),whichClusCrysBelongsTo_.end(),seedHit.id(),PairSortByFirst<DetId,int>());
213 
214  if(result.first!=result.second) protoClusters_[result.first->second].removeHit(seedHit);
215  protoClusters_[clusNr].addSeed();
216 
217  }
218  }
219  }
220 
221 
222 
223  for(size_t clusNr=0;clusNr<protoClusters_.size();clusNr++){
224  const ProtoBasicCluster& protoCluster= protoClusters_[clusNr];
225  Point position;
226  position = posCalculator_.Calculate_Location(protoCluster.hits(), hits, geometry_p, geometryES_p);
227  clusters_v.push_back(reco::BasicCluster(protoCluster.energy(), position, reco::CaloID(detector_), protoCluster.hits(),
228  reco::CaloCluster::multi5x5, protoCluster.seed().id()));
229  }
230 
231  protoClusters_.clear();
232  whichClusCrysBelongsTo_.clear();
233 }
234 
237  const CaloSubdetectorGeometry *geometryES,
239  bool seedOutside)
240 {
241 
242  double energy = 0;
243  //double chi2 = 0;
244  reco::CaloID caloID;
245  Point position;
246  position = posCalculator_.Calculate_Location(current_v, hits, geometry, geometryES);
247 
248  std::vector<std::pair<DetId, float> >::iterator it;
249  for (it = current_v.begin(); it != current_v.end(); it++)
250  {
251  EcalRecHitCollection::const_iterator itt = hits->find( (*it).first );
252  EcalRecHit hit_p = *itt;
253  energy += hit_p.energy();
254  //chi2 += 0;
255  if ( (*it).first.subdetId() == EcalBarrel ) {
257  } else {
259  }
260 
261  }
262  //chi2 /= energy;
263 
264  LogTrace("EcalClusters") << "******** NEW CLUSTER ********";
265  LogTrace("EcalClusters") << "No. of crystals = " << current_v.size();
266  LogTrace("EcalClusters") << " Energy = " << energy ;
267  LogTrace("EcalClusters") << " Phi = " << position.phi();
268  LogTrace("EcalClusters") << " Eta " << position.eta();
269  LogTrace("EcalClusters") << "*****************************";
270 
271 
272  // to be a valid cluster the cluster energy
273  // must be at least the seed energy
274  double seedEnergy = seedIt->energy();
275  if ((seedOutside && energy>=0) || (!seedOutside && energy >= seedEnergy))
276  {
277  if(reassignSeedCrysToClusterItSeeds_){ //if we're not doing this, we dont need this info so lets not bother filling it
278  for(size_t hitNr=0;hitNr<current_v.size();hitNr++) whichClusCrysBelongsTo_.push_back(std::pair<DetId,int>(current_v[hitNr].first,protoClusters_.size()));
279  }
280  protoClusters_.push_back(ProtoBasicCluster(energy,*seedIt,current_v));
281 
282  // clusters_v.push_back(reco::BasicCluster(energy, position, reco::CaloID(detector_), current_v, reco::CaloCluster::multi5x5, seedIt->id()));
283 
284  // if no valid cluster was built,
285  // then free up these crystals to be used in the next...
286  } else {
287  std::vector<std::pair<DetId, float> >::iterator iter;
288  for (iter = current_v.begin(); iter != current_v.end(); iter++)
289  {
290  used_s.erase(iter->first);
291  } //for(iter)
292  } //else
293 
294 }
295 
297  const EcalRecHitCollection *hits)
298 {
299 
300  bool maxima = true;
302  EcalRecHitCollection::const_iterator seedHit = hits->find(navigator.pos());
303  double seedEnergy = seedHit->energy();
304 
305  std::vector<DetId> swissCrossVec;
306  swissCrossVec.clear();
307 
308  swissCrossVec.push_back(navigator.west());
309  navigator.home();
310  swissCrossVec.push_back(navigator.east());
311  navigator.home();
312  swissCrossVec.push_back(navigator.north());
313  navigator.home();
314  swissCrossVec.push_back(navigator.south());
315  navigator.home();
316 
317  std::vector<DetId>::const_iterator detItr;
318  for (unsigned int i = 0; i < swissCrossVec.size(); ++i)
319  {
320 
321  // look for this hit
322  thisHit = recHits_->find(swissCrossVec[i]);
323 
324  // continue if this hit was not found
325  if ((swissCrossVec[i] == DetId(0)) || thisHit == recHits_->end()) continue;
326 
327  // the recHit has to be skipped in the local maximum search if it was found
328  // in the map of channels to be excluded
329  uint32_t rhFlag = thisHit->recoFlag();
330  std::vector<int>::const_iterator vit = std::find(v_chstatus_.begin(), v_chstatus_.end(), rhFlag);
331  if (vit != v_chstatus_.end()) continue;
332 
333  // if this crystal has more energy than the seed then we do
334  // not have a local maxima
335  if (thisHit->energy() > seedEnergy)
336  {
337  maxima = false;
338  break;
339  }
340  }
341 
342  return maxima;
343 
344 }
345 
347  const EcalRecHitCollection *hits,
349 {
350 
351  DetId thisDet;
352  std::set<DetId>::iterator setItr;
353 
354  // now add the 5x5 taking care to mark the edges
355  // as able to seed and where overlapping in the central
356  // region with crystals that were previously able to seed
357  // change their status so they are not able to seed
358  //std::cout << std::endl;
359  for (int dx = -2; dx < 3; ++dx)
360  {
361  for (int dy = -2; dy < 3; ++ dy)
362  {
363 
364  // navigate in free steps forming
365  // a full 5x5
366  thisDet = navigator.offsetBy(dx, dy);
367  navigator.home();
368 
369  // add the current crystal
370  //std::cout << "adding " << dx << ", " << dy << std::endl;
371  addCrystal(thisDet);
372 
373  // now consider if we are in an edge (outer 16)
374  // or central (inner 9) region
375  if ((abs(dx) > 1) || (abs(dy) > 1))
376  {
377  // this is an "edge" so should be allowed to seed
378  // provided it is not already used
379  //std::cout << " setting can seed" << std::endl;
380  canSeed_s.insert(thisDet);
381  } // end if "edge"
382  else
383  {
384  // or else we are in the central 3x3
385  // and must remove any of these crystals from the canSeed set
386  setItr = canSeed_s.find(thisDet);
387  if (setItr != canSeed_s.end())
388  {
389  //std::cout << " unsetting can seed" << std::endl;
390  canSeed_s.erase(setItr);
391  }
392  } // end if "centre"
393 
394 
395  } // end loop on dy
396 
397  } // end loop on dx
398 
399  //std::cout << "*** " << std::endl;
400  //std::cout << " current_v contains " << current_v.size() << std::endl;
401  //std::cout << "*** " << std::endl;
402 }
403 
404 
406 {
407 
408  EcalRecHitCollection::const_iterator thisIt = recHits_->find(det);
409  if ((thisIt != recHits_->end()) && (thisIt->id() != DetId(0)))
410  {
411  if ((used_s.find(thisIt->id()) == used_s.end()))
412  {
413  //std::cout << " ... this is a good crystal and will be added" << std::endl;
414  current_v.push_back( std::pair<DetId, float>(det, 1.) ); // by default hit energy fractions are set at 1.
415  used_s.insert(det);
416  }
417  }
418 
419 }
420 
422 {
423  std::vector<std::pair<DetId,float> >::iterator hitPos;
424  for(hitPos=hits_.begin();hitPos<hits_.end();hitPos++){
425  if(hitToRM.id()==hitPos->first) break;
426  }
427  if(hitPos!=hits_.end()){
428  hits_.erase(hitPos);
429  energy_-=hitToRM.energy();
430  return true;
431  }return false;
432 }
433 
434 //now the tricky thing is to make sure we insert it in the right place in the hit vector
435 //order is from -2,-2, -2,-1 to 2,2 with seed being 0,0 (eta,phi)
437 {
438  typedef std::vector<std::pair<DetId,float> >::iterator It;
439  std::pair<It,It> hitPos;
440 
441  if(seed_.id().subdetId()==EcalBarrel){
442  hitPos = std::equal_range(hits_.begin(),hits_.end(),seed_.id(),PairSortByFirst<DetId,float,EBDetIdSorter>());
443  }else{
444  hitPos = std::equal_range(hits_.begin(),hits_.end(),seed_.id(),PairSortByFirst<DetId,float,EEDetIdSorter>());
445  }
446 
447  if(hitPos.first==hitPos.second){//it doesnt already exist in the vec, add it
448  hits_.insert(hitPos.first,std::pair<DetId,float>(seed_.id(),1.));
449  energy_+=seed_.energy();
450  containsSeed_=true;
451 
452  return true;
453  }else return false;
454 
455 }
456 
458 {
459  for(size_t hitNr=0;hitNr<hits_.size();hitNr++){
460  if(seed_.id()==hits_[hitNr].first) return true;
461  }
462  return false;
463 }
int ix() const
Definition: EEDetId.h:76
void makeCluster(const EcalRecHitCollection *hits, const CaloSubdetectorGeometry *geometry_p, const CaloSubdetectorGeometry *geometryES_p, const EcalRecHitCollection::const_iterator &seedIt, bool seedOutside)
std::vector< reco::BasicCluster > makeClusters(const EcalRecHitCollection *hits, const CaloSubdetectorGeometry *geometry, const CaloSubdetectorTopology *topology_p, const CaloSubdetectorGeometry *geometryES_p, reco::CaloID::Detectors detector, bool regional=false, const std::vector< EcalEtaPhiRegion > &regions=std::vector< EcalEtaPhiRegion >())
std::vector< EcalRecHit >::const_iterator const_iterator
const std::vector< std::pair< DetId, float > > & hits() const
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:20
int iphi() const
get the crystal iphi
Definition: EBDetId.h:53
void mainSearch(const EcalRecHitCollection *hits, const CaloSubdetectorGeometry *geometry_p, const CaloSubdetectorTopology *topology_p, const CaloSubdetectorGeometry *geometryES_p)
T offsetBy(int deltaX, int deltaY) const
Free movement of arbitray steps.
Definition: CaloNavigator.h:80
void prepareCluster(CaloNavigator< DetId > &navigator, const EcalRecHitCollection *hits, const CaloSubdetectorGeometry *geometry)
bool operator()(const std::pair< T1, T2 > &lhs, const T1 &rhs)
T west() const
move the navigator west
Definition: CaloNavigator.h:59
void addCrystal(const DetId &det)
int zside() const
Definition: EEDetId.h:70
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
float energy() const
Definition: EcalRecHit.h:68
bool operator()(const EEDetId &lhs, const EEDetId &rhs)
int iy() const
Definition: EEDetId.h:82
T south() const
move the navigator south
Definition: CaloNavigator.h:45
int ieta() const
get the crystal ieta
Definition: EBDetId.h:51
T pos() const
get the current position
Definition: CaloNavigator.h:32
bool operator()(const T1 &lhs, const T1 &rhs)
#define LogTrace(id)
const_iterator end() const
bool removeHit(const EcalRecHit &hitToRM)
T east() const
move the navigator east
Definition: CaloNavigator.h:52
void home() const
move the navigator back to the starting point
Definition: DetId.h:18
DetId id() const
get the id
Definition: EcalRecHit.h:77
bool operator()(const EBDetId &lhs, const EBDetId &rhs)
bool operator()(const std::pair< T1, T2 > &lhs, const std::pair< T1, T2 > &rhs)
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.
bool checkMaxima(CaloNavigator< DetId > &navigator, const EcalRecHitCollection *hits)
bool operator()(const T1 &lhs, const std::pair< T1, T2 > &rhs)
math::XYZPoint Point
point in the space
iterator find(key_type k)
static int position[264][3]
Definition: ReadPGInfo.cc:509
T north() const
move the navigator north
Definition: CaloNavigator.h:38
#define ET
const_iterator begin() const