CMS 3D CMS Logo

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

#include <HybridClusterAlgo.h>

Public Member Functions

 HybridClusterAlgo ()
 
 HybridClusterAlgo (double eb_str, int step, double ethres, double eseed, double xi, bool useEtForXi, double ewing, const std::vector< int > &v_chstatus, const PositionCalc &posCalculator, bool dynamicEThres=false, double eThresA=0, double eThresB=0.1, const std::vector< int > &severityToExclude=std::vector< int >().operator=(std::vector< int >(999)), bool excludeFromCluster=false)
 
void mainSearch (const EcalRecHitCollection *hits, const CaloSubdetectorGeometry *geometry)
 
void makeClusters (const EcalRecHitCollection *, const CaloSubdetectorGeometry *geometry, reco::BasicClusterCollection &basicClusters, const EcalSeverityLevelAlgo *sevLv, bool regional=false, const std::vector< RectangularEtaPhiRegion > &regions=std::vector< RectangularEtaPhiRegion >())
 
double makeDomino (EcalBarrelNavigatorHT &navigator, std::vector< EcalRecHit > &cells)
 
reco::SuperClusterCollection makeSuperClusters (const reco::CaloClusterPtrVector &)
 
void setDynamicPhiRoad (const edm::ParameterSet &bremRecoveryPset)
 
 ~HybridClusterAlgo ()
 

Private Types

typedef math::XYZPoint Point
 

Private Member Functions

double e2Et (EcalBarrelNavigatorHT &navigator, const EcalRecHitCollection *hits, const CaloSubdetectorGeometry *geometry)
 
double et25 (EcalBarrelNavigatorHT &navigator, const EcalRecHitCollection *hits, const CaloSubdetectorGeometry *geometry)
 

Private Attributes

std::map< int, std::vector< reco::BasicCluster > > clustered_
 
bool dynamicEThres_
 
bool dynamicPhiRoad_
 
double eb_st
 
double Eseed
 
double eThres_
 
double eThresA_
 
double eThresB_
 
double Ewing
 
std::set< DetIdexcludedCrys_
 
bool excludeFromCluster_
 
BremRecoveryPhiRoadAlgophiRoadAlgo_
 
int phiSteps_
 
PositionCalc posCalculator_
 
const EcalRecHitCollectionrecHits_
 
std::vector< reco::BasicClusterseedClus_
 
std::vector< EcalRecHitseeds
 
float severityRecHitThreshold_
 
float severitySpikeThreshold_
 
EcalBarrelHardcodedTopologytopo_
 
std::set< DetIduseddetids
 
bool UseEtForXi
 
std::vector< int > v_chstatus_
 
std::vector< int > v_severitylevel_
 
double Xi
 

Detailed Description

Definition at line 26 of file HybridClusterAlgo.h.

Member Typedef Documentation

◆ Point

Definition at line 29 of file HybridClusterAlgo.h.

Constructor & Destructor Documentation

◆ HybridClusterAlgo() [1/2]

HybridClusterAlgo::HybridClusterAlgo ( )
inline

Definition at line 113 of file HybridClusterAlgo.h.

113 {}

◆ HybridClusterAlgo() [2/2]

HybridClusterAlgo::HybridClusterAlgo ( double  eb_str,
int  step,
double  ethres,
double  eseed,
double  xi,
bool  useEtForXi,
double  ewing,
const std::vector< int > &  v_chstatus,
const PositionCalc posCalculator,
bool  dynamicEThres = false,
double  eThresA = 0,
double  eThresB = 0.1,
const std::vector< int > &  severityToExclude = std::vector<int>().operator=(std::vector<int>(999)),
bool  excludeFromCluster = false 
)

Definition at line 15 of file HybridClusterAlgo.cc.

References dynamicEThres_, dynamicPhiRoad_, Eseed, eThresA_, eThresB_, LogTrace, posCalculator_, jetsAK4_CHS_cff::sort, topo_, UseEtForXi, v_chstatus_, and Xi.

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),
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 }
EcalBarrelHardcodedTopology * topo_
#define LogTrace(id)
PositionCalc posCalculator_
std::vector< int > v_chstatus_
std::vector< int > v_severitylevel_
step
Definition: StallMonitor.cc:98

◆ ~HybridClusterAlgo()

HybridClusterAlgo::~HybridClusterAlgo ( )
inline

Definition at line 136 of file HybridClusterAlgo.h.

References dynamicPhiRoad_, phiRoadAlgo_, and topo_.

136  {
137  if (dynamicPhiRoad_)
138  delete phiRoadAlgo_;
139  delete topo_;
140  // delete SCShape_;
141  }
EcalBarrelHardcodedTopology * topo_
BremRecoveryPhiRoadAlgo * phiRoadAlgo_

Member Function Documentation

◆ e2Et()

double HybridClusterAlgo::e2Et ( EcalBarrelNavigatorHT navigator,
const EcalRecHitCollection hits,
const CaloSubdetectorGeometry geometry 
)
private

Definition at line 640 of file HybridClusterAlgo.cc.

References PositionCalc::Calculate_Location(), PVValHelper::dx, PVValHelper::dy, edm::SortedCollection< T, SORT >::end(), edm::SortedCollection< T, SORT >::find(), hfClusterShapes_cfi::hits, HLT_2022v15_cff::navigator, posCalculator_, and recHits_.

Referenced by mainSearch().

642  {
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
664  return 1 / cosh(pos.eta());
665 }
std::vector< EcalRecHit >::const_iterator const_iterator
PositionCalc posCalculator_
math::XYZPoint Calculate_Location(const HitsAndFractions &iDetIds, const edm::SortedCollection< HitType > *iRecHits, const CaloSubdetectorGeometry *iSubGeom, const CaloSubdetectorGeometry *iESGeom=nullptr)
Definition: PositionCalc.h:65
const_iterator end() const
Definition: DetId.h:17
iterator find(key_type k)
Structure Point Contains parameters of Gaussian fits to DMRs.
const EcalRecHitCollection * recHits_

◆ et25()

double HybridClusterAlgo::et25 ( EcalBarrelNavigatorHT navigator,
const EcalRecHitCollection hits,
const CaloSubdetectorGeometry geometry 
)
private

Definition at line 608 of file HybridClusterAlgo.cc.

References PositionCalc::Calculate_Location(), PVValHelper::dx, PVValHelper::dy, edm::SortedCollection< T, SORT >::end(), CastorDataFrameFilter_impl::energySum(), EgHLTOffHistBins_cfi::et, edm::SortedCollection< T, SORT >::find(), hfClusterShapes_cfi::hits, HLT_2022v15_cff::navigator, posCalculator_, and recHits_.

Referenced by mainSearch().

610  {
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;
636  double et = energySum / cosh(pos.eta());
637  return et;
638 }
std::vector< EcalRecHit >::const_iterator const_iterator
PositionCalc posCalculator_
math::XYZPoint Calculate_Location(const HitsAndFractions &iDetIds, const edm::SortedCollection< HitType > *iRecHits, const CaloSubdetectorGeometry *iSubGeom, const CaloSubdetectorGeometry *iESGeom=nullptr)
Definition: PositionCalc.h:65
const_iterator end() const
Definition: DetId.h:17
iterator find(key_type k)
Structure Point Contains parameters of Gaussian fits to DMRs.
const EcalRecHitCollection * recHits_
double energySum(const DataFrame &df, int fs, int ls)

◆ mainSearch()

void HybridClusterAlgo::mainSearch ( const EcalRecHitCollection hits,
const CaloSubdetectorGeometry geometry 
)

Definition at line 173 of file HybridClusterAlgo.cc.

References BremRecoveryPhiRoadAlgo::barrelPhiRoad(), PositionCalc::Calculate_Location(), clustered_, reco::CaloID::DET_ECAL_BARREL, dynamicEThres_, dynamicPhiRoad_, e2Et(), Eseed, et25(), eThres_, eThresA_, eThresB_, hfClusterShapes_cfi::hits, reco::CaloCluster::hybrid, mps_fire::i, createfilelist::int, dqmiolumiharvest::j, dqmdumpme::k, LogTrace, makeDomino(), HLT_2022v15_cff::navigator, nhits, DetId::null(), phiRoadAlgo_, phiSteps_, posCalculator_, FastTrackerRecHitMaskProducer_cfi::recHits, seedClus_, seeds, mathSSE::sqrt(), groupFilesInBlocks::temp, topo_, useddetids, UseEtForXi, and Xi.

Referenced by makeClusters().

173  {
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
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(
427  // and also add to the vector of clusters that will be used in constructing
428  // the supercluster
429  thisseedClusters.push_back(reco::BasicCluster(
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
EcalBarrelHardcodedTopology * topo_
std::map< int, std::vector< reco::BasicCluster > > clustered_
double e2Et(EcalBarrelNavigatorHT &navigator, const EcalRecHitCollection *hits, const CaloSubdetectorGeometry *geometry)
#define LogTrace(id)
PositionCalc posCalculator_
constexpr bool null() const
is this a null id ?
Definition: DetId.h:59
std::vector< EcalRecHit > seeds
T sqrt(T t)
Definition: SSEVec.h:19
BremRecoveryPhiRoadAlgo * phiRoadAlgo_
math::XYZPoint Calculate_Location(const HitsAndFractions &iDetIds, const edm::SortedCollection< HitType > *iRecHits, const CaloSubdetectorGeometry *iSubGeom, const CaloSubdetectorGeometry *iESGeom=nullptr)
Definition: PositionCalc.h:65
std::set< DetId > useddetids
double et25(EcalBarrelNavigatorHT &navigator, const EcalRecHitCollection *hits, const CaloSubdetectorGeometry *geometry)
Definition: DetId.h:17
double makeDomino(EcalBarrelNavigatorHT &navigator, std::vector< EcalRecHit > &cells)
Structure Point Contains parameters of Gaussian fits to DMRs.
std::vector< reco::BasicCluster > seedClus_

◆ makeClusters()

void HybridClusterAlgo::makeClusters ( const EcalRecHitCollection recColl,
const CaloSubdetectorGeometry geometry,
reco::BasicClusterCollection basicClusters,
const EcalSeverityLevelAlgo sevLv,
bool  regional = false,
const std::vector< RectangularEtaPhiRegion > &  regions = std::vector<RectangularEtaPhiRegion>() 
)

Definition at line 60 of file HybridClusterAlgo.cc.

References edm::SortedCollection< T, SORT >::begin(), clustered_, eb_st, edm::SortedCollection< T, SORT >::end(), ET, excludedCrys_, excludeFromCluster_, spr::find(), createfilelist::int, isClusterEtLess(), dqmiolumiharvest::j, EcalRecHit::kGood, LogTrace, mainSearch(), position, recHits_, HLT_2022v15_cff::region, seedClus_, seeds, EcalSeverityLevelAlgo::severityLevel(), jetsAK4_CHS_cff::sort, useddetids, v_chstatus_, v_severitylevel_, x, and y.

Referenced by HybridClusterProducer::produce().

65  {
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 }
void mainSearch(const EcalRecHitCollection *hits, const CaloSubdetectorGeometry *geometry)
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
#define LogTrace(id)
EcalSeverityLevel::SeverityLevel severityLevel(const DetId &id) const
Evaluate status from id use channelStatus from DB.
std::vector< EcalRecHit > seeds
std::vector< int > v_chstatus_
const_iterator begin() const
bool isClusterEtLess(const reco::CaloCluster &x, const reco::CaloCluster &y)
Definition: ClusterEtLess.h:7
std::set< DetId > useddetids
std::set< DetId > excludedCrys_
const_iterator end() const
std::vector< int > v_severitylevel_
std::vector< BasicCluster > BasicClusterCollection
collection of BasicCluster objects
static int position[264][3]
Definition: ReadPGInfo.cc:289
const EcalRecHitCollection * recHits_
#define ET
std::vector< reco::BasicCluster > seedClus_

◆ makeDomino()

double HybridClusterAlgo::makeDomino ( EcalBarrelNavigatorHT navigator,
std::vector< EcalRecHit > &  cells 
)

Definition at line 525 of file HybridClusterAlgo.cc.

References cells, edm::SortedCollection< T, SORT >::end(), EcalRecHit::energy(), Ewing, excludedCrys_, edm::SortedCollection< T, SORT >::find(), hcaltpdigi_cfi::ieta1, hcaltpdigi_cfi::ieta2, hcaltpdigi_cfi::ieta3, hcaltpdigi_cfi::ieta4, HLT_2022v15_cff::navigator, recHits_, and useddetids.

Referenced by mainSearch().

525  {
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 }
std::vector< EcalRecHit >::const_iterator const_iterator
std::set< DetId > useddetids
std::set< DetId > excludedCrys_
const_iterator end() const
caConstants::TupleMultiplicity const CAHitNtupletGeneratorKernelsGPU::HitToTuple const cms::cuda::AtomicPairCounter GPUCACell const *__restrict__ cells
Definition: DetId.h:17
iterator find(key_type k)
const EcalRecHitCollection * recHits_
float energy() const
Definition: EcalRecHit.h:68

◆ makeSuperClusters()

reco::SuperClusterCollection HybridClusterAlgo::makeSuperClusters ( const reco::CaloClusterPtrVector clustersCollection)

Definition at line 450 of file HybridClusterAlgo.cc.

References clustered_, reco::CaloCluster::energy(), mps_fire::i, createfilelist::int, isClusterEtLess(), dqmiolumiharvest::j, LogTrace, reco::CaloCluster::position(), RecoTauValidation_cfi::posX, RecoTauValidation_cfi::posY, edm::PtrVector< T >::push_back(), fileCollector::seed, seedClus_, edm::PtrVectorBase::size(), and jetsAK4_CHS_cff::sort.

Referenced by HybridClusterProducer::produce().

450  {
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 }
const math::XYZPoint & position() const
cluster centroid position
Definition: CaloCluster.h:154
size_type size() const
Size of the RefVector.
Definition: PtrVectorBase.h:75
void push_back(Ptr< T > const &iPtr)
Definition: PtrVector.h:149
std::map< int, std::vector< reco::BasicCluster > > clustered_
#define LogTrace(id)
std::vector< SuperCluster > SuperClusterCollection
collection of SuperCluser objectr
bool isClusterEtLess(const reco::CaloCluster &x, const reco::CaloCluster &y)
Definition: ClusterEtLess.h:7
double energy() const
cluster energy
Definition: CaloCluster.h:149
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
std::vector< reco::BasicCluster > seedClus_

◆ setDynamicPhiRoad()

void HybridClusterAlgo::setDynamicPhiRoad ( const edm::ParameterSet bremRecoveryPset)
inline

Member Data Documentation

◆ clustered_

std::map<int, std::vector<reco::BasicCluster> > HybridClusterAlgo::clustered_
private

Definition at line 95 of file HybridClusterAlgo.h.

Referenced by mainSearch(), makeClusters(), and makeSuperClusters().

◆ dynamicEThres_

bool HybridClusterAlgo::dynamicEThres_
private

Definition at line 71 of file HybridClusterAlgo.h.

Referenced by HybridClusterAlgo(), and mainSearch().

◆ dynamicPhiRoad_

bool HybridClusterAlgo::dynamicPhiRoad_
private

◆ eb_st

double HybridClusterAlgo::eb_st
private

Definition at line 32 of file HybridClusterAlgo.h.

Referenced by makeClusters().

◆ Eseed

double HybridClusterAlgo::Eseed
private

Definition at line 56 of file HybridClusterAlgo.h.

Referenced by HybridClusterAlgo(), and mainSearch().

◆ eThres_

double HybridClusterAlgo::eThres_
private

Definition at line 51 of file HybridClusterAlgo.h.

Referenced by mainSearch().

◆ eThresA_

double HybridClusterAlgo::eThresA_
private

Definition at line 52 of file HybridClusterAlgo.h.

Referenced by HybridClusterAlgo(), and mainSearch().

◆ eThresB_

double HybridClusterAlgo::eThresB_
private

Definition at line 53 of file HybridClusterAlgo.h.

Referenced by HybridClusterAlgo(), and mainSearch().

◆ Ewing

double HybridClusterAlgo::Ewing
private

Definition at line 65 of file HybridClusterAlgo.h.

Referenced by makeDomino().

◆ excludedCrys_

std::set<DetId> HybridClusterAlgo::excludedCrys_
private

Definition at line 109 of file HybridClusterAlgo.h.

Referenced by makeClusters(), and makeDomino().

◆ excludeFromCluster_

bool HybridClusterAlgo::excludeFromCluster_
private

Definition at line 108 of file HybridClusterAlgo.h.

Referenced by makeClusters().

◆ phiRoadAlgo_

BremRecoveryPhiRoadAlgo* HybridClusterAlgo::phiRoadAlgo_
private

Definition at line 48 of file HybridClusterAlgo.h.

Referenced by mainSearch(), setDynamicPhiRoad(), and ~HybridClusterAlgo().

◆ phiSteps_

int HybridClusterAlgo::phiSteps_
private

Definition at line 37 of file HybridClusterAlgo.h.

Referenced by mainSearch().

◆ posCalculator_

PositionCalc HybridClusterAlgo::posCalculator_
private

Definition at line 98 of file HybridClusterAlgo.h.

Referenced by e2Et(), et25(), HybridClusterAlgo(), and mainSearch().

◆ recHits_

const EcalRecHitCollection* HybridClusterAlgo::recHits_
private

Definition at line 78 of file HybridClusterAlgo.h.

Referenced by e2Et(), et25(), makeClusters(), and makeDomino().

◆ seedClus_

std::vector<reco::BasicCluster> HybridClusterAlgo::seedClus_
private

Definition at line 92 of file HybridClusterAlgo.h.

Referenced by mainSearch(), makeClusters(), and makeSuperClusters().

◆ seeds

std::vector<EcalRecHit> HybridClusterAlgo::seeds
private

Definition at line 89 of file HybridClusterAlgo.h.

Referenced by mainSearch(), and makeClusters().

◆ severityRecHitThreshold_

float HybridClusterAlgo::severityRecHitThreshold_
private

Definition at line 105 of file HybridClusterAlgo.h.

◆ severitySpikeThreshold_

float HybridClusterAlgo::severitySpikeThreshold_
private

Definition at line 106 of file HybridClusterAlgo.h.

◆ topo_

EcalBarrelHardcodedTopology* HybridClusterAlgo::topo_
private

Definition at line 81 of file HybridClusterAlgo.h.

Referenced by HybridClusterAlgo(), mainSearch(), and ~HybridClusterAlgo().

◆ useddetids

std::set<DetId> HybridClusterAlgo::useddetids
private

Definition at line 86 of file HybridClusterAlgo.h.

Referenced by mainSearch(), makeClusters(), and makeDomino().

◆ UseEtForXi

bool HybridClusterAlgo::UseEtForXi
private

Definition at line 62 of file HybridClusterAlgo.h.

Referenced by HybridClusterAlgo(), and mainSearch().

◆ v_chstatus_

std::vector<int> HybridClusterAlgo::v_chstatus_
private

Definition at line 101 of file HybridClusterAlgo.h.

Referenced by HybridClusterAlgo(), and makeClusters().

◆ v_severitylevel_

std::vector<int> HybridClusterAlgo::v_severitylevel_
private

Definition at line 104 of file HybridClusterAlgo.h.

Referenced by makeClusters().

◆ Xi

double HybridClusterAlgo::Xi
private

Definition at line 59 of file HybridClusterAlgo.h.

Referenced by HybridClusterAlgo(), and mainSearch().