CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
HybridClusterAlgo.cc
Go to the documentation of this file.
8 
9 #include <iostream>
10 #include <map>
11 #include <vector>
12 #include <set>
13 
14 //The real constructor
16  int step,
17  double ethres,
18  double eseed,
19  double xi,
20  bool useEtForXi,
21  double ewing,
22  const std::vector<int>& v_chstatus,
23  const PositionCalc& posCalculator,
24  bool dynamicEThres,
25  double eThresA,
26  double eThresB,
27  const std::vector<int>& severityToExclude,
28  bool excludeFromCluster)
29  :
30 
31  eb_st(eb_str),
32  phiSteps_(step),
33  eThres_(ethres),
34  eThresA_(eThresA),
35  eThresB_(eThresB),
36  Eseed(eseed),
37  Xi(xi),
38  UseEtForXi(useEtForXi),
39  Ewing(ewing),
40  dynamicEThres_(dynamicEThres),
41  v_chstatus_(v_chstatus),
42  v_severitylevel_(severityToExclude),
43  //severityRecHitThreshold_(severityRecHitThreshold),
44  //severitySpikeThreshold_(severitySpikeThreshold),
45  excludeFromCluster_(excludeFromCluster)
46 
47 {
48  dynamicPhiRoad_ = false;
49  LogTrace("EcalClusters") << "dynamicEThres: " << dynamicEThres_ << " : A,B " << eThresA_ << ", " << eThresB_;
50  LogTrace("EcalClusters") << "Eseed: " << Eseed << " , Xi: " << Xi << ", UseEtForXi: " << UseEtForXi;
51 
52  //if (dynamicPhiRoad_) phiRoadAlgo_ = new BremRecoveryPhiRoadAlgo(bremRecoveryPset);
53  posCalculator_ = posCalculator;
55 
56  std::sort(v_chstatus_.begin(), v_chstatus_.end());
57 }
58 
59 // Return a vector of clusters from a collection of EcalRecHits:
62  reco::BasicClusterCollection& basicClusters,
63  const EcalSeverityLevelAlgo* sevLv,
64  bool regional,
65  const std::vector<RectangularEtaPhiRegion>& regions) {
66  //clear vector of seeds
67  seeds.clear();
68  //clear flagged crystals
69  excludedCrys_.clear();
70  //clear map of supercluster/basiccluster association
71  clustered_.clear();
72  //clear set of used detids
73  useddetids.clear();
74  //clear vector of seed clusters
75  seedClus_.clear();
76  //Pass in a pointer to the collection.
77  recHits_ = recColl;
78 
79  LogTrace("EcalClusters") << "Cleared vectors, starting clusterization...";
80  LogTrace("EcalClusters") << " Purple monkey aardvark.";
81  //
82 
83  int nregions = 0;
84  if (regional)
85  nregions = regions.size();
86 
87  if (!regional || nregions) {
89 
90  for (it = recHits_->begin(); it != recHits_->end(); it++) {
91  //Make the vector of seeds that we're going to use.
92  //One of the few places position is used, needed for ET calculation.
93  auto this_cell = geometry->getGeometry(it->id());
94  const GlobalPoint& position = this_cell->getPosition();
95 
96  // Require that RecHit is within clustering region in case
97  // of regional reconstruction
98  bool withinRegion = false;
99  if (regional) {
100  std::vector<RectangularEtaPhiRegion>::const_iterator region;
101  for (region = regions.begin(); region != regions.end(); region++) {
102  if (region->inRegion(this_cell->etaPos(), this_cell->phiPos())) {
103  withinRegion = true;
104  break;
105  }
106  }
107  }
108 
109  if (!regional || withinRegion) {
110  //Must pass seed threshold
111  float ET = it->energy() * position.basicVector().unit().perp();
112 
113  LogTrace("EcalClusters") << "Seed crystal: ";
114 
115  if (ET > eb_st) {
116  // avoid seeding for anomalous channels
117  if (!it->checkFlag(EcalRecHit::kGood)) { // if rechit is good, no need for further checks
118  if (it->checkFlags(v_chstatus_)) {
120  excludedCrys_.insert(it->id());
121  continue; // the recHit has to be excluded from seeding
122  }
123  }
124 
125  int severityFlag = sevLv->severityLevel(it->id(), *recHits_);
126  std::vector<int>::const_iterator sit =
127  std::find(v_severitylevel_.begin(), v_severitylevel_.end(), severityFlag);
128  LogTrace("EcalClusters") << "found flag: " << severityFlag;
129 
130  if (sit != v_severitylevel_.end()) {
132  excludedCrys_.insert(it->id());
133  continue;
134  }
135  seeds.push_back(*it);
136 
137  LogTrace("EcalClusters") << "Seed ET: " << ET;
138  LogTrace("EcalClsuters") << "Seed E: " << it->energy();
139  }
140  }
141  }
142  }
143 
144  //Yay sorting.
145  LogTrace("EcalClusters") << "Built vector of seeds, about to sort them...";
146 
147  //Needs three argument sort with seed comparison operator
148  sort(seeds.begin(), seeds.end(), [](auto const& x, auto const& y) { return x.energy() > y.energy(); });
149 
150  LogTrace("EcalClusters") << "done";
151 
152  //Now to do the work.
153  LogTrace("EcalClusters") << "About to call mainSearch...";
154  mainSearch(recColl, geometry);
155  LogTrace("EcalClusters") << "done";
156 
157  //Hand the basicclusters back to the producer. It has to
158  //put them in the event. Then we can make superclusters.
159  std::map<int, reco::BasicClusterCollection>::iterator bic;
160  for (bic = clustered_.begin(); bic != clustered_.end(); bic++) {
161  reco::BasicClusterCollection bl = bic->second;
162  for (int j = 0; j < int(bl.size()); ++j) {
163  basicClusters.push_back(bl[j]);
164  }
165  }
166 
167  //Yay more sorting.
168  sort(basicClusters.rbegin(), basicClusters.rend(), isClusterEtLess);
169  //Done!
170  LogTrace("EcalClusters") << "returning to producer. ";
171 }
172 
174  LogTrace("EcalClusters") << "HybridClusterAlgo Algorithm - looking for clusters";
175  LogTrace("EcalClusters") << "Found the following clusters:";
176 
177  // Loop over seeds:
178  std::vector<EcalRecHit>::iterator it;
179  int clustercounter = 0;
180 
181  for (it = seeds.begin(); it != seeds.end(); it++) {
182  std::vector<reco::BasicCluster> thisseedClusters;
183  DetId itID = it->id();
184 
185  // make sure the current seed has not been used/will not be used in the future:
186  std::set<DetId>::iterator seed_in_rechits_it = useddetids.find(itID);
187 
188  if (seed_in_rechits_it != useddetids.end())
189  continue;
190  //If this seed is already used, then don't use it again.
191 
192  // output some info on the hit:
193  LogTrace("EcalClusters") << "*****************************************************"
194  << "Seed of energy E = " << it->energy() << " @ " << EBDetId(itID)
195  << "*****************************************************";
196 
197  //Make a navigator, and set it to the seed cell.
199 
200  //Now use the navigator to start building dominoes.
201 
202  //Walking positive in phi:
203  std::vector<double> dominoEnergyPhiPlus; //Here I will store the results of the domino sums
204  std::vector<std::vector<EcalRecHit> > dominoCellsPhiPlus; //These are the actual EcalRecHit for dominos.
205 
206  //Walking negative in phi
207  std::vector<double> dominoEnergyPhiMinus; //Here I will store the results of the domino sums
208  std::vector<std::vector<EcalRecHit> > dominoCellsPhiMinus; //These are the actual EcalRecHit for dominos.
209 
210  //The two sets together.
211  std::vector<double> dominoEnergy; //Here I will store the results of the domino sums
212  std::vector<std::vector<EcalRecHit> > dominoCells; //These are the actual EcalRecHit for dominos.
213 
214  //First, the domino about the seed:
215  std::vector<EcalRecHit> initialdomino;
216  double e_init = makeDomino(navigator, initialdomino);
217  if (e_init < Eseed)
218  continue;
219 
220  LogTrace("EcalClusters") << "Made initial domino";
221 
222  //
223  // compute the phi road length
224  double phiSteps;
225  double et5x5 = et25(navigator, hits, geometry);
226  if (dynamicPhiRoad_ && e_init > 0) {
227  phiSteps = phiRoadAlgo_->barrelPhiRoad(et5x5);
228  navigator.home();
229  } else
230  phiSteps = phiSteps_;
231 
232  //Positive phi steps.
233  for (int i = 0; i < phiSteps; ++i) {
234  //remember, this always increments the current position of the navigator.
235  DetId centerD = navigator.north();
236  if (centerD.null())
237  continue;
238  LogTrace("EcalClusters") << "Step ++" << i << " @ " << EBDetId(centerD);
239  EcalBarrelNavigatorHT dominoNav(centerD, topo_);
240 
241  //Go get the new domino.
242  std::vector<EcalRecHit> dcells;
243  double etemp = makeDomino(dominoNav, dcells);
244 
245  //save this information
246  dominoEnergyPhiPlus.push_back(etemp);
247  dominoCellsPhiPlus.push_back(dcells);
248  }
249  LogTrace("EcalClusters") << "Got positive dominos";
250  //return to initial position
251  navigator.home();
252 
253  //Negative phi steps.
254  for (int i = 0; i < phiSteps; ++i) {
255  //remember, this always decrements the current position of the navigator.
256  DetId centerD = navigator.south();
257  if (centerD.null())
258  continue;
259 
260  LogTrace("EcalClusters") << "Step --" << i << " @ " << EBDetId(centerD);
261  EcalBarrelNavigatorHT dominoNav(centerD, topo_);
262 
263  //Go get the new domino.
264  std::vector<EcalRecHit> dcells;
265  double etemp = makeDomino(dominoNav, dcells);
266 
267  //save this information
268  dominoEnergyPhiMinus.push_back(etemp);
269  dominoCellsPhiMinus.push_back(dcells);
270  }
271 
272  LogTrace("EcalClusters") << "Got negative dominos: ";
273 
274  //Assemble this information:
275  for (int i = int(dominoEnergyPhiMinus.size()) - 1; i >= 0; --i) {
276  dominoEnergy.push_back(dominoEnergyPhiMinus[i]);
277  dominoCells.push_back(dominoCellsPhiMinus[i]);
278  }
279  dominoEnergy.push_back(e_init);
280  dominoCells.push_back(initialdomino);
281  for (int i = 0; i < int(dominoEnergyPhiPlus.size()); ++i) {
282  dominoEnergy.push_back(dominoEnergyPhiPlus[i]);
283  dominoCells.push_back(dominoCellsPhiPlus[i]);
284  }
285 
286  //Ok, now I have all I need in order to go ahead and make clusters.
287  LogTrace("EcalClusters") << "Dumping domino energies: ";
288 
289  //for (int i=0;i<int(dominoEnergy.size());++i){
290  // LogTrace("EcalClusters") << "Domino: " << i << " E: " << dominoEnergy[i] ;
291  // }
292  //Identify the peaks in this set of dominos:
293  //Peak==a domino whose energy is greater than the two adjacent dominos.
294  //thus a peak in the local sense.
295  double etToE(1.);
296  if (!UseEtForXi) {
297  etToE = 1. / e2Et(navigator, hits, geometry);
298  }
299  std::vector<int> PeakIndex;
300  for (int i = 1; i < int(dominoEnergy.size()) - 1; ++i) {
301  if (dominoEnergy[i] > dominoEnergy[i - 1] && dominoEnergy[i] >= dominoEnergy[i + 1] &&
302  dominoEnergy[i] > sqrt(Eseed * Eseed + Xi * Xi * et5x5 * et5x5 * etToE *
303  etToE)) { // Eseed increases ~proportiaonally to et5x5
304  PeakIndex.push_back(i);
305  }
306  }
307 
308  LogTrace("EcalClusters") << "Found: " << PeakIndex.size() << " peaks.";
309 
310  //Order these peaks by energy:
311  for (int i = 0; i < int(PeakIndex.size()); ++i) {
312  for (int j = 0; j < int(PeakIndex.size()) - 1; ++j) {
313  if (dominoEnergy[PeakIndex[j]] < dominoEnergy[PeakIndex[j + 1]]) {
314  int ihold = PeakIndex[j + 1];
315  PeakIndex[j + 1] = PeakIndex[j];
316  PeakIndex[j] = ihold;
317  }
318  }
319  }
320 
321  std::vector<int> OwnerShip;
322  std::vector<double> LumpEnergy;
323  OwnerShip.reserve(int(dominoEnergy.size()));
324  for (int i = 0; i < int(dominoEnergy.size()); ++i)
325  OwnerShip.push_back(-1);
326 
327  //Loop over peaks.
328  double eThres = eThres_;
329  double e5x5 = 0.0;
330  for (int i = 0; i < int(PeakIndex.size()); ++i) {
331  int idxPeak = PeakIndex[i];
332  OwnerShip[idxPeak] = i;
333  double lump = dominoEnergy[idxPeak];
334 
335  // compute eThres for this peak
336  // if set to dynamic (otherwise uncanged from
337  // fixed setting
338  if (dynamicEThres_) {
339  //std::cout << "i : " << i << " idxPeak " << idxPeak << std::endl;
340  //std::cout << " the dominoEnergy.size() = " << dominoEnergy.size() << std::endl;
341  // compute e5x5 for this seed crystal
342  //std::cout << "idxPeak, phiSteps " << idxPeak << ", " << phiSteps << std::endl;
343  e5x5 = lump;
344  //std::cout << "lump " << e5x5 << std::endl;
345  if ((idxPeak + 1) < (int)dominoEnergy.size())
346  e5x5 += dominoEnergy[idxPeak + 1];
347  //std::cout << "+1 " << e5x5 << std::endl;
348  if ((idxPeak + 2) < (int)dominoEnergy.size())
349  e5x5 += dominoEnergy[idxPeak + 2];
350  //std::cout << "+2 " << e5x5 << std::endl;
351  if ((idxPeak - 1) > 0)
352  e5x5 += dominoEnergy[idxPeak - 1];
353  //std::cout << "-1 " << e5x5 << std::endl;
354  if ((idxPeak - 2) > 0)
355  e5x5 += dominoEnergy[idxPeak - 2];
356  //std::cout << "-2 " << e5x5 << std::endl;
357  // compute eThres
358  eThres = (eThresA_ * e5x5) + eThresB_;
359  //std::cout << eThres << std::endl;
360  //std::cout << std::endl;
361  }
362 
363  //Loop over adjacent dominos at higher phi
364  for (int j = idxPeak + 1; j < int(dominoEnergy.size()); ++j) {
365  if (OwnerShip[j] == -1 && dominoEnergy[j] > eThres && dominoEnergy[j] <= dominoEnergy[j - 1]) {
366  OwnerShip[j] = i;
367  lump += dominoEnergy[j];
368  } else {
369  break;
370  }
371  }
372  //loop over adjacent dominos at lower phi. Sum up energy of lumps.
373  for (int j = idxPeak - 1; j >= 0; --j) {
374  if (OwnerShip[j] == -1 && dominoEnergy[j] > eThres && dominoEnergy[j] <= dominoEnergy[j + 1]) {
375  OwnerShip[j] = i;
376  lump += dominoEnergy[j];
377  } else {
378  break;
379  }
380  }
381  LumpEnergy.push_back(lump);
382  }
383 
384  //Make the basic clusters:
385  for (int i = 0; i < int(PeakIndex.size()); ++i) {
386  bool HasSeedCrystal = false;
387  //One cluster for each peak.
388  std::vector<EcalRecHit> recHits;
389  std::vector<std::pair<DetId, float> > dets;
390  int nhits = 0;
391  for (int j = 0; j < int(dominoEnergy.size()); ++j) {
392  if (OwnerShip[j] == i) {
393  std::vector<EcalRecHit> temp = dominoCells[j];
394  for (int k = 0; k < int(temp.size()); ++k) {
395  dets.push_back(std::pair<DetId, float>(temp[k].id(), 1.)); // by default energy fractions are 1
396  if (temp[k].id() == itID)
397  HasSeedCrystal = true;
398  recHits.push_back(temp[k]);
399  nhits++;
400  }
401  }
402  }
403  LogTrace("EcalClusters") << "Adding a cluster with: " << nhits;
404  LogTrace("EcalClusters") << "total E: " << LumpEnergy[i];
405  LogTrace("EcalClusters") << "total dets: " << dets.size();
406 
407  //Get Calorimeter position
408  Point pos = posCalculator_.Calculate_Location(dets, hits, geometry);
409 
410  //double totChi2=0;
411  //double totE=0;
412  std::vector<std::pair<DetId, float> > usedHits;
413  for (int blarg = 0; blarg < int(recHits.size()); ++blarg) {
414  //totChi2 +=0;
415  //totE+=recHits[blarg].energy();
416  usedHits.push_back(std::make_pair<DetId, float>(recHits[blarg].id(), 1.0));
417  useddetids.insert(recHits[blarg].id());
418  }
419 
420  //if (totE>0)
421  //totChi2/=totE;
422 
423  if (HasSeedCrystal) {
424  // note that this "basiccluster" has the seed crystal of the hyrbid, so record it
425  seedClus_.push_back(reco::BasicCluster(
426  LumpEnergy[i], pos, reco::CaloID(reco::CaloID::DET_ECAL_BARREL), usedHits, reco::CaloCluster::hybrid, itID));
427  // and also add to the vector of clusters that will be used in constructing
428  // the supercluster
429  thisseedClusters.push_back(reco::BasicCluster(
430  LumpEnergy[i], pos, reco::CaloID(reco::CaloID::DET_ECAL_BARREL), usedHits, reco::CaloCluster::hybrid, itID));
431  } else {
432  // note that if this "basiccluster" is not the one that seeded the hybrid,
433  // the seed crystal is unset in this entry in the vector of clusters that will
434  // be used in constructing the super cluster
435  thisseedClusters.push_back(reco::BasicCluster(
437  }
438  }
439 
440  // Make association so that superclusters can be made later.
441  // but only if some BasicClusters have been found...
442  if (!thisseedClusters.empty()) {
443  clustered_.insert(std::make_pair(clustercounter, thisseedClusters));
444  clustercounter++;
445  }
446  } //Seed loop
447 
448 } // end of mainSearch
449 
451  //Here's what we'll return.
453 
454  //Here's our map iterator that gives us the appropriate association.
455  std::map<int, reco::BasicClusterCollection>::iterator mapit;
456  for (mapit = clustered_.begin(); mapit != clustered_.end(); mapit++) {
458  reco::CaloClusterPtr seed; //This is not really a seed, but I need to tell SuperCluster something.
459  //So I choose the highest energy basiccluster in the SuperCluster.
460 
461  std::vector<reco::BasicCluster> thiscoll = mapit->second; //This is the set of BasicClusters in this
462  //SuperCluster
463 
464  double ClusterE = 0; //Sum of cluster energies for supercluster.
465  //Holders for position of this supercluster.
466  double posX = 0;
467  double posY = 0;
468  double posZ = 0;
469 
470  //Loop over this set of basic clusters, find their references, and add them to the
471  //supercluster. This could be somehow more efficient.
472 
473  for (int i = 0; i < int(thiscoll.size()); ++i) {
474  reco::BasicCluster thisclus = thiscoll[i]; //The Cluster in question.
475  for (int j = 0; j < int(clustersCollection.size()); ++j) {
476  //Find the appropriate cluster from the list of references
477  reco::BasicCluster cluster_p = *clustersCollection[j];
478  if (thisclus == cluster_p) { //Comparison based on energy right now.
479  thissc.push_back(clustersCollection[j]);
480  bool isSeed = false;
481  for (int qu = 0; qu < int(seedClus_.size()); ++qu) {
482  if (cluster_p == seedClus_[qu])
483  isSeed = true;
484  }
485  if (isSeed)
486  seed = clustersCollection[j];
487 
488  ClusterE += cluster_p.energy();
489  posX += cluster_p.energy() * cluster_p.position().X();
490  posY += cluster_p.energy() * cluster_p.position().Y();
491  posZ += cluster_p.energy() * cluster_p.position().Z();
492  }
493  } //End loop over finding references.
494  } //End loop over clusters.
495 
496  posX /= ClusterE;
497  posY /= ClusterE;
498  posZ /= ClusterE;
499 
500  /* //This part is moved to EgammaSCEnergyCorrectionAlgo
501  double preshowerE = 0.;
502  double phiWidth = 0.;
503  double etaWidth = 0.;
504  //Calculate phiWidth & etaWidth for SuperClusters
505  reco::SuperCluster suCl(ClusterE, math::XYZPoint(posX, posY, posZ), seed, thissc, preshowerE, phiWidth, etaWidth);
506  SCShape_->Calculate_Covariances(suCl);
507  phiWidth = SCShape_->phiWidth();
508  etaWidth = SCShape_->etaWidth();
509  //Assign phiWidth & etaWidth to SuperCluster as data members
510  suCl.setPhiWidth(phiWidth);
511  suCl.setEtaWidth(etaWidth);
512  */
513 
514  reco::SuperCluster suCl(ClusterE, math::XYZPoint(posX, posY, posZ), seed, thissc);
515  SCcoll.push_back(suCl);
516 
517  LogTrace("EcalClusters") << "Super cluster sum: " << ClusterE;
518  LogTrace("EcalClusters") << "Made supercluster with energy E: " << suCl.energy();
519 
520  } //end loop over map
521  sort(SCcoll.rbegin(), SCcoll.rend(), isClusterEtLess);
522  return SCcoll;
523 }
524 
526  //At the beginning of this function, the navigator starts at the middle of the domino,
527  //and that's where EcalBarrelNavigator::home() should send it.
528  //Walk one crystal in eta to either side of the initial point. Sum the three cell energies.
529  //If the resultant energy is larger than Ewing, then go an additional cell to either side.
530  //Returns: Total domino energy. Also, stores the cells used to create domino in the vector.
531  cells.clear();
532  double Etot = 0;
533 
534  //Ready? Get the starting cell.
535  DetId center = navigator.pos();
537 
538  if (center_it != recHits_->end()) {
539  EcalRecHit SeedHit = *center_it;
540  if (useddetids.find(center) == useddetids.end() && excludedCrys_.find(center) == excludedCrys_.end()) {
541  Etot += SeedHit.energy();
542  cells.push_back(SeedHit);
543  }
544  }
545  //One step upwards in Ieta:
546  DetId ieta1 = navigator.west();
548  if (eta1_it != recHits_->end()) {
549  EcalRecHit UpEta = *eta1_it;
550  if (useddetids.find(ieta1) == useddetids.end() && excludedCrys_.find(ieta1) == excludedCrys_.end()) {
551  Etot += UpEta.energy();
552  cells.push_back(UpEta);
553  }
554  }
555 
556  //Go back to the middle.
557  navigator.home();
558 
559  //One step downwards in Ieta:
560  DetId ieta2 = navigator.east();
562  if (eta2_it != recHits_->end()) {
563  EcalRecHit DownEta = *eta2_it;
564  if (useddetids.find(ieta2) == useddetids.end() && excludedCrys_.find(ieta2) == excludedCrys_.end()) {
565  Etot += DownEta.energy();
566  cells.push_back(DownEta);
567  }
568  }
569 
570  //Now check the energy. If smaller than Ewing, then we're done. If greater than Ewing, we have to
571  //add two additional cells, the 'wings'
572  if (Etot < Ewing) {
573  navigator.home(); //Needed even here!!
574  return Etot; //Done! Not adding 'wings'.
575  }
576 
577  //Add the extra 'wing' cells. Remember, we haven't sent the navigator home,
578  //we're still on the DownEta cell.
579  if (ieta2 != DetId(0)) {
580  DetId ieta3 = navigator.east(); //Take another step downward.
582  if (eta3_it != recHits_->end()) {
583  EcalRecHit DownEta2 = *eta3_it;
584  if (useddetids.find(ieta3) == useddetids.end() && excludedCrys_.find(ieta3) == excludedCrys_.end()) {
585  Etot += DownEta2.energy();
586  cells.push_back(DownEta2);
587  }
588  }
589  }
590  //Now send the navigator home.
591  navigator.home();
592  navigator.west(); //Now you're on eta1_it
593  if (ieta1 != DetId(0)) {
594  DetId ieta4 = navigator.west(); //Take another step upward.
596  if (eta4_it != recHits_->end()) {
597  EcalRecHit UpEta2 = *eta4_it;
598  if (useddetids.find(ieta4) == useddetids.end() && excludedCrys_.find(ieta4) == excludedCrys_.end()) {
599  Etot += UpEta2.energy();
600  cells.push_back(UpEta2);
601  }
602  }
603  }
604  navigator.home();
605  return Etot;
606 }
607 
609  const EcalRecHitCollection* hits,
611  DetId thisDet;
612  std::vector<std::pair<DetId, float> > dets;
613  dets.clear();
615  double energySum = 0.0;
616 
617  for (int dx = -2; dx < 3; ++dx) {
618  for (int dy = -2; dy < 3; ++dy) {
619  //std::cout << "dx, dy " << dx << ", " << dy << std::endl;
620  thisDet = navigator.offsetBy(dx, dy);
621  navigator.home();
622 
623  if (thisDet != DetId(0)) {
624  hit = recHits_->find(thisDet);
625  if (hit != recHits_->end()) {
626  dets.push_back(std::pair<DetId, float>(thisDet, 1.)); // by default hit energy fraction is set to 1
627  energySum += hit->energy();
628  }
629  }
630  }
631  }
632 
633  // convert it to ET
634  //std::cout << "dets.size(), energySum: " << dets.size() << ", " << energySum << std::endl;
635  Point pos = posCalculator_.Calculate_Location(dets, hits, geometry);
636  double et = energySum / cosh(pos.eta());
637  return et;
638 }
639 
641  const EcalRecHitCollection* hits,
643  DetId thisDet;
644  std::vector<std::pair<DetId, float> > dets;
645  dets.clear();
647 
648  for (int dx = -2; dx < 3; ++dx) {
649  for (int dy = -2; dy < 3; ++dy) {
650  thisDet = navigator.offsetBy(dx, dy);
651  navigator.home();
652 
653  if (thisDet != DetId(0)) {
654  hit = recHits_->find(thisDet);
655  if (hit != recHits_->end()) {
656  dets.push_back(std::pair<DetId, float>(thisDet, 1.)); // by default hit energy fraction is set to 1
657  }
658  }
659  }
660  }
661 
662  // compute coefficient to turn E into Et
663  Point pos = posCalculator_.Calculate_Location(dets, hits, geometry);
664  return 1 / cosh(pos.eta());
665 }
const math::XYZPoint & position() const
cluster centroid position
Definition: CaloCluster.h:154
size_type size() const
Size of the RefVector.
Definition: PtrVectorBase.h:75
EcalSeverityLevel::SeverityLevel severityLevel(const DetId &id) const
Evaluate status from id use channelStatus from DB.
constexpr bool null() const
is this a null id ?
Definition: DetId.h:59
void push_back(Ptr< T > const &iPtr)
Definition: PtrVector.h:149
EcalBarrelHardcodedTopology * topo_
void mainSearch(const EcalRecHitCollection *hits, const CaloSubdetectorGeometry *geometry)
Basic3DVector unit() const
std::map< int, std::vector< reco::BasicCluster > > clustered_
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
double e2Et(EcalBarrelNavigatorHT &navigator, const EcalRecHitCollection *hits, const CaloSubdetectorGeometry *geometry)
#define LogTrace(id)
PositionCalc posCalculator_
T offsetBy(int deltaX, int deltaY) const
Free movement of arbitray steps.
Definition: CaloNavigator.h:66
std::vector< EcalRecHit > seeds
std::vector< SuperCluster > SuperClusterCollection
collection of SuperCluser objectr
T west() const
move the navigator west
Definition: CaloNavigator.h:48
T sqrt(T t)
Definition: SSEVec.h:19
reco::SuperClusterCollection makeSuperClusters(const reco::CaloClusterPtrVector &)
BremRecoveryPhiRoadAlgo * phiRoadAlgo_
float energy() const
Definition: EcalRecHit.h:68
math::XYZPoint Calculate_Location(const HitsAndFractions &iDetIds, const edm::SortedCollection< HitType > *iRecHits, const CaloSubdetectorGeometry *iSubGeom, const CaloSubdetectorGeometry *iESGeom=nullptr)
Definition: PositionCalc.h:65
double energy() const
cluster energy
Definition: CaloCluster.h:149
std::vector< int > v_chstatus_
T south() const
move the navigator south
Definition: CaloNavigator.h:36
bool isClusterEtLess(const reco::CaloCluster &x, const reco::CaloCluster &y)
Definition: ClusterEtLess.h:7
std::set< DetId > useddetids
T pos() const
get the current position
Definition: CaloNavigator.h:24
double et25(EcalBarrelNavigatorHT &navigator, const EcalRecHitCollection *hits, const CaloSubdetectorGeometry *geometry)
T perp() const
Magnitude of transverse component.
std::set< DetId > excludedCrys_
const_iterator end() const
std::vector< int > v_severitylevel_
caConstants::TupleMultiplicity const CAHitNtupletGeneratorKernelsGPU::HitToTuple const cms::cuda::AtomicPairCounter GPUCACell const *__restrict__ cells
T east() const
move the navigator east
Definition: CaloNavigator.h:42
void home() const
move the navigator back to the starting point
Definition: CaloNavigator.h:96
Definition: DetId.h:17
double makeDomino(EcalBarrelNavigatorHT &navigator, std::vector< EcalRecHit > &cells)
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.
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
std::vector< BasicCluster > BasicClusterCollection
collection of BasicCluster objects
iterator find(key_type k)
Structure Point Contains parameters of Gaussian fits to DMRs.
Definition: DMRtrends.cc:57
static int position[264][3]
Definition: ReadPGInfo.cc:289
T north() const
move the navigator north
Definition: CaloNavigator.h:30
const EcalRecHitCollection * recHits_
#define ET
step
Definition: StallMonitor.cc:94
double energySum(const DataFrame &df, int fs, int ls)
void makeClusters(const EcalRecHitCollection *, const CaloSubdetectorGeometry *geometry, reco::BasicClusterCollection &basicClusters, const EcalSeverityLevelAlgo *sevLv, bool regional=false, const std::vector< RectangularEtaPhiRegion > &regions=std::vector< RectangularEtaPhiRegion >())
const BasicVectorType & basicVector() const
Definition: PV3DBase.h:53
std::vector< reco::BasicCluster > seedClus_
const_iterator begin() const