CMS 3D CMS Logo

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