22 const std::vector<int>& v_chstatus,
27 const std::vector<int>& severityToExclude,
28 bool excludeFromCluster)
38 UseEtForXi(useEtForXi),
40 dynamicEThres_(dynamicEThres),
41 v_chstatus_(v_chstatus),
42 v_severitylevel_(severityToExclude),
45 excludeFromCluster_(excludeFromCluster)
65 const std::vector<RectangularEtaPhiRegion>&
regions) {
79 LogTrace(
"EcalClusters") <<
"Cleared vectors, starting clusterization...";
80 LogTrace(
"EcalClusters") <<
" Purple monkey aardvark.";
85 nregions = regions.size();
87 if (!regional || nregions) {
98 bool withinRegion =
false;
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())) {
109 if (!regional || withinRegion) {
113 LogTrace(
"EcalClusters") <<
"Seed crystal: ";
126 std::vector<int>::const_iterator sit =
128 LogTrace(
"EcalClusters") <<
"found flag: " << severityFlag;
135 seeds.push_back(*it);
137 LogTrace(
"EcalClusters") <<
"Seed ET: " <<
ET;
138 LogTrace(
"EcalClsuters") <<
"Seed E: " << it->energy();
145 LogTrace(
"EcalClusters") <<
"Built vector of seeds, about to sort them...";
148 sort(
seeds.begin(),
seeds.end(), [](
auto const&
x,
auto const&
y) {
return x.energy() >
y.energy(); });
153 LogTrace(
"EcalClusters") <<
"About to call mainSearch...";
159 std::map<int, reco::BasicClusterCollection>::iterator bic;
162 for (
int j = 0;
j < int(bl.size()); ++
j) {
163 basicClusters.push_back(bl[
j]);
170 LogTrace(
"EcalClusters") <<
"returning to producer. ";
174 LogTrace(
"EcalClusters") <<
"HybridClusterAlgo Algorithm - looking for clusters";
175 LogTrace(
"EcalClusters") <<
"Found the following clusters:";
178 std::vector<EcalRecHit>::iterator it;
179 int clustercounter = 0;
181 for (it =
seeds.begin(); it !=
seeds.end(); it++) {
182 std::vector<reco::BasicCluster> thisseedClusters;
183 DetId itID = it->id();
186 std::set<DetId>::iterator seed_in_rechits_it =
useddetids.find(itID);
193 LogTrace(
"EcalClusters") <<
"*****************************************************"
194 <<
"Seed of energy E = " << it->energy() <<
" @ " <<
EBDetId(itID)
195 <<
"*****************************************************";
203 std::vector<double> dominoEnergyPhiPlus;
204 std::vector<std::vector<EcalRecHit> > dominoCellsPhiPlus;
207 std::vector<double> dominoEnergyPhiMinus;
208 std::vector<std::vector<EcalRecHit> > dominoCellsPhiMinus;
211 std::vector<double> dominoEnergy;
212 std::vector<std::vector<EcalRecHit> > dominoCells;
215 std::vector<EcalRecHit> initialdomino;
216 double e_init =
makeDomino(navigator, initialdomino);
220 LogTrace(
"EcalClusters") <<
"Made initial domino";
225 double et5x5 =
et25(navigator, hits, geometry);
233 for (
int i = 0;
i < phiSteps; ++
i) {
242 std::vector<EcalRecHit> dcells;
246 dominoEnergyPhiPlus.push_back(etemp);
247 dominoCellsPhiPlus.push_back(dcells);
249 LogTrace(
"EcalClusters") <<
"Got positive dominos";
254 for (
int i = 0;
i < phiSteps; ++
i) {
264 std::vector<EcalRecHit> dcells;
268 dominoEnergyPhiMinus.push_back(etemp);
269 dominoCellsPhiMinus.push_back(dcells);
272 LogTrace(
"EcalClusters") <<
"Got negative dominos: ";
275 for (
int i =
int(dominoEnergyPhiMinus.size()) - 1;
i >= 0; --
i) {
276 dominoEnergy.push_back(dominoEnergyPhiMinus[
i]);
277 dominoCells.push_back(dominoCellsPhiMinus[i]);
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]);
287 LogTrace(
"EcalClusters") <<
"Dumping domino energies: ";
297 etToE = 1. /
e2Et(navigator, hits, geometry);
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] &&
304 PeakIndex.push_back(
i);
308 LogTrace(
"EcalClusters") <<
"Found: " << PeakIndex.size() <<
" peaks.";
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;
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);
330 for (
int i = 0;
i < int(PeakIndex.size()); ++
i) {
331 int idxPeak = PeakIndex[
i];
332 OwnerShip[idxPeak] =
i;
333 double lump = dominoEnergy[idxPeak];
345 if ((idxPeak + 1) < (int)dominoEnergy.size())
346 e5x5 += dominoEnergy[idxPeak + 1];
348 if ((idxPeak + 2) < (int)dominoEnergy.size())
349 e5x5 += dominoEnergy[idxPeak + 2];
351 if ((idxPeak - 1) > 0)
352 e5x5 += dominoEnergy[idxPeak - 1];
354 if ((idxPeak - 2) > 0)
355 e5x5 += dominoEnergy[idxPeak - 2];
364 for (
int j = idxPeak + 1;
j < int(dominoEnergy.size()); ++
j) {
365 if (OwnerShip[
j] == -1 && dominoEnergy[
j] > eThres && dominoEnergy[
j] <= dominoEnergy[
j - 1]) {
367 lump += dominoEnergy[
j];
373 for (
int j = idxPeak - 1;
j >= 0; --
j) {
374 if (OwnerShip[
j] == -1 && dominoEnergy[
j] > eThres && dominoEnergy[
j] <= dominoEnergy[
j + 1]) {
376 lump += dominoEnergy[
j];
381 LumpEnergy.push_back(lump);
385 for (
int i = 0;
i < int(PeakIndex.size()); ++
i) {
386 bool HasSeedCrystal =
false;
388 std::vector<EcalRecHit>
recHits;
389 std::vector<std::pair<DetId, float> > dets;
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.));
396 if (temp[
k].
id() == itID)
397 HasSeedCrystal =
true;
398 recHits.push_back(temp[
k]);
403 LogTrace(
"EcalClusters") <<
"Adding a cluster with: " <<
nhits;
404 LogTrace(
"EcalClusters") <<
"total E: " << LumpEnergy[
i];
405 LogTrace(
"EcalClusters") <<
"total dets: " << dets.size();
412 std::vector<std::pair<DetId, float> > usedHits;
413 for (
int blarg = 0; blarg < int(recHits.size()); ++blarg) {
416 usedHits.push_back(std::make_pair<DetId, float>(recHits[blarg].
id(), 1.0));
423 if (HasSeedCrystal) {
442 if (!thisseedClusters.empty()) {
443 clustered_.insert(std::make_pair(clustercounter, thisseedClusters));
455 std::map<int, reco::BasicClusterCollection>::iterator mapit;
461 std::vector<reco::BasicCluster> thiscoll = mapit->second;
473 for (
int i = 0;
i < int(thiscoll.size()); ++
i) {
475 for (
int j = 0;
j < int(clustersCollection.
size()); ++
j) {
478 if (thisclus == cluster_p) {
481 for (
int qu = 0; qu < int(
seedClus_.size()); ++qu) {
486 seed = clustersCollection[
j];
488 ClusterE += cluster_p.
energy();
515 SCcoll.push_back(suCl);
517 LogTrace(
"EcalClusters") <<
"Super cluster sum: " << ClusterE;
518 LogTrace(
"EcalClusters") <<
"Made supercluster with energy E: " << suCl.
energy();
542 cells.push_back(SeedHit);
552 cells.push_back(UpEta);
566 cells.push_back(DownEta);
579 if (ieta2 !=
DetId(0)) {
585 Etot += DownEta2.
energy();
586 cells.push_back(DownEta2);
593 if (ieta1 !=
DetId(0)) {
600 cells.push_back(UpEta2);
612 std::vector<std::pair<DetId, float> > dets;
617 for (
int dx = -2;
dx < 3; ++
dx) {
618 for (
int dy = -2;
dy < 3; ++
dy) {
623 if (thisDet !=
DetId(0)) {
626 dets.push_back(std::pair<DetId, float>(thisDet, 1.));
627 energySum += hit->energy();
636 double et = energySum / cosh(pos.eta());
644 std::vector<std::pair<DetId, float> > dets;
648 for (
int dx = -2;
dx < 3; ++
dx) {
649 for (
int dy = -2;
dy < 3; ++
dy) {
653 if (thisDet !=
DetId(0)) {
656 dets.push_back(std::pair<DetId, float>(thisDet, 1.));
664 return 1 / cosh(pos.eta());
const math::XYZPoint & position() const
cluster centroid position
size_type size() const
Size of the RefVector.
int barrelPhiRoad(double et)
EcalSeverityLevel::SeverityLevel severityLevel(const DetId &id) const
Evaluate status from id use channelStatus from DB.
constexpr bool null() const
is this a null id ?
void push_back(Ptr< T > const &iPtr)
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)
double e2Et(EcalBarrelNavigatorHT &navigator, const EcalRecHitCollection *hits, const CaloSubdetectorGeometry *geometry)
PositionCalc posCalculator_
T offsetBy(int deltaX, int deltaY) const
Free movement of arbitray steps.
std::vector< EcalRecHit > seeds
std::vector< SuperCluster > SuperClusterCollection
collection of SuperCluser objectr
T west() const
move the navigator west
reco::SuperClusterCollection makeSuperClusters(const reco::CaloClusterPtrVector &)
BremRecoveryPhiRoadAlgo * phiRoadAlgo_
math::XYZPoint Calculate_Location(const HitsAndFractions &iDetIds, const edm::SortedCollection< HitType > *iRecHits, const CaloSubdetectorGeometry *iSubGeom, const CaloSubdetectorGeometry *iESGeom=nullptr)
double energy() const
cluster energy
std::vector< int > v_chstatus_
T south() const
move the navigator south
bool isClusterEtLess(const reco::CaloCluster &x, const reco::CaloCluster &y)
std::set< DetId > useddetids
T pos() const
get the current position
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
void home() const
move the navigator back to the starting point
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
std::vector< BasicCluster > BasicClusterCollection
collection of BasicCluster objects
iterator find(key_type k)
Structure Point Contains parameters of Gaussian fits to DMRs.
static int position[264][3]
T north() const
move the navigator north
const EcalRecHitCollection * recHits_
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 > ®ions=std::vector< RectangularEtaPhiRegion >())
const BasicVectorType & basicVector() const
std::vector< reco::BasicCluster > seedClus_
const_iterator begin() const