23 const std::vector<int>& v_chstatus,
28 const std::vector<int>& severityToExclude,
29 bool excludeFromCluster
32 eb_st(eb_str), phiSteps_(step),
33 eThres_(ethres), eThresA_(eThresA), eThresB_(eThresB),
34 Eseed(eseed), Xi(xi), UseEtForXi(useEtForXi), Ewing(ewing),
35 dynamicEThres_(dynamicEThres),
36 v_chstatus_(v_chstatus), v_severitylevel_(severityToExclude),
39 excludeFromCluster_(excludeFromCluster)
63 const std::vector<RectangularEtaPhiRegion>&
regions 79 LogTrace(
"EcalClusters") <<
"Cleared vectors, starting clusterization..." ;
80 LogTrace(
"EcalClusters") <<
" Purple monkey aardvark." ;
84 if(regional) nregions=regions.size();
86 if(!regional || nregions) {
100 bool withinRegion =
false;
102 std::vector<RectangularEtaPhiRegion>::const_iterator region;
103 for (region=regions.begin(); region!=regions.end(); region++) {
104 if (region->inRegion(this_cell->etaPos(),this_cell->phiPos())) {
111 if (!regional || withinRegion) {
116 LogTrace(
"EcalClusters") <<
"Seed crystal: " ;
132 LogTrace(
"EcalClusters") <<
"found flag: " << severityFlag ;
140 seeds.push_back(*it);
142 LogTrace(
"EcalClusters") <<
"Seed ET: " <<
ET ;
143 LogTrace(
"EcalClsuters") <<
"Seed E: " << it->energy() ;
152 LogTrace(
"EcalClusters") <<
"Built vector of seeds, about to sort them..." ;
155 sort(
seeds.begin(),
seeds.end(), [](
auto const&
x,
auto const&
y) {
return x.energy() >
y.energy() ; });
157 LogTrace(
"EcalClusters") <<
"done" ;
160 LogTrace(
"EcalClusters") <<
"About to call mainSearch...";
162 LogTrace(
"EcalClusters") <<
"done" ;
166 std::map<int, reco::BasicClusterCollection>::iterator bic;
169 for (
int j=0;j<
int(bl.size());++j){
170 basicClusters.push_back(bl[j]);
177 LogTrace(
"EcalClusters") <<
"returning to producer. ";
184 LogTrace(
"EcalClusters") <<
"HybridClusterAlgo Algorithm - looking for clusters" ;
185 LogTrace(
"EcalClusters") <<
"Found the following clusters:" ;
188 std::vector<EcalRecHit>::iterator it;
189 int clustercounter=0;
191 for (it =
seeds.begin(); it !=
seeds.end(); it++){
192 std::vector <reco::BasicCluster> thisseedClusters;
193 DetId itID = it->id();
196 std::set<DetId>::iterator seed_in_rechits_it =
useddetids.find(itID);
198 if (seed_in_rechits_it !=
useddetids.end())
continue;
202 LogTrace(
"EcalClusters") <<
"*****************************************************" <<
"Seed of energy E = " << it->energy() <<
" @ " <<
EBDetId(itID) <<
"*****************************************************" ;
211 std::vector <double> dominoEnergyPhiPlus;
212 std::vector <std::vector <EcalRecHit> > dominoCellsPhiPlus;
215 std::vector <double> dominoEnergyPhiMinus;
216 std::vector <std::vector <EcalRecHit> > dominoCellsPhiMinus;
219 std::vector <double> dominoEnergy;
220 std::vector <std::vector <EcalRecHit> > dominoCells;
223 std::vector <EcalRecHit> initialdomino;
224 double e_init =
makeDomino(navigator, initialdomino);
225 if (e_init <
Eseed)
continue;
227 LogTrace(
"EcalClusters") <<
"Made initial domino" ;
232 double et5x5 =
et25(navigator, hits, geometry);
240 for (
int i=0;
i<phiSteps;++
i){
249 std::vector <EcalRecHit> dcells;
253 dominoEnergyPhiPlus.push_back(etemp);
254 dominoCellsPhiPlus.push_back(dcells);
256 LogTrace(
"EcalClusters") <<
"Got positive dominos" ;
261 for (
int i=0;
i<phiSteps;++
i){
271 std::vector <EcalRecHit> dcells;
275 dominoEnergyPhiMinus.push_back(etemp);
276 dominoCellsPhiMinus.push_back(dcells);
279 LogTrace(
"EcalClusters") <<
"Got negative dominos: " ;
282 for (
int i=
int(dominoEnergyPhiMinus.size())-1;
i >= 0;--
i){
283 dominoEnergy.push_back(dominoEnergyPhiMinus[
i]);
284 dominoCells.push_back(dominoCellsPhiMinus[i]);
286 dominoEnergy.push_back(e_init);
287 dominoCells.push_back(initialdomino);
288 for (
int i=0;
i<
int(dominoEnergyPhiPlus.size());++
i){
289 dominoEnergy.push_back(dominoEnergyPhiPlus[
i]);
290 dominoCells.push_back(dominoCellsPhiPlus[i]);
294 LogTrace(
"EcalClusters") <<
"Dumping domino energies: " ;
305 etToE = 1./
e2Et(navigator, hits, geometry);
307 std::vector <int> PeakIndex;
308 for (
int i=1;
i<
int(dominoEnergy.size())-1;++
i){
309 if (dominoEnergy[
i] > dominoEnergy[
i-1]
310 && dominoEnergy[
i] >= dominoEnergy[
i+1]
313 PeakIndex.push_back(
i);
317 LogTrace(
"EcalClusters") <<
"Found: " << PeakIndex.size() <<
" peaks." ;
320 for (
int i=0;
i<
int(PeakIndex.size());++
i){
321 for (
int j=0;j<
int(PeakIndex.size())-1;++j){
322 if (dominoEnergy[PeakIndex[j]] < dominoEnergy[PeakIndex[j+1]]){
323 int ihold = PeakIndex[j+1];
324 PeakIndex[j+1] = PeakIndex[j];
325 PeakIndex[j] = ihold;
330 std::vector<int> OwnerShip;
331 std::vector<double> LumpEnergy;
332 for (
int i=0;
i<
int(dominoEnergy.size());++
i) OwnerShip.push_back(-1);
337 for (
int i = 0;
i <
int(PeakIndex.size()); ++
i)
340 int idxPeak = PeakIndex[
i];
341 OwnerShip[idxPeak] =
i;
342 double lump = dominoEnergy[idxPeak];
355 if ((idxPeak + 1) < (
int)dominoEnergy.size()) e5x5 += dominoEnergy[idxPeak + 1];
357 if ((idxPeak + 2) < (
int)dominoEnergy.size()) e5x5 += dominoEnergy[idxPeak + 2];
359 if ((idxPeak - 1) > 0) e5x5 += dominoEnergy[idxPeak - 1];
361 if ((idxPeak - 2) > 0) e5x5 += dominoEnergy[idxPeak - 2];
370 for (
int j=idxPeak+1;j<
int(dominoEnergy.size());++j){
371 if (OwnerShip[j]==-1 &&
372 dominoEnergy[j] > eThres
373 && dominoEnergy[j] <= dominoEnergy[j-1]){
375 lump+=dominoEnergy[j];
382 for (
int j=idxPeak-1;j>=0;--j){
383 if (OwnerShip[j]==-1 &&
384 dominoEnergy[j] > eThres
385 && dominoEnergy[j] <= dominoEnergy[j+1]){
387 lump+=dominoEnergy[j];
393 LumpEnergy.push_back(lump);
397 for (
int i=0;
i<
int(PeakIndex.size());++
i){
398 bool HasSeedCrystal =
false;
400 std::vector<EcalRecHit> recHits;
401 std::vector< std::pair<DetId, float> > dets;
403 for (
int j=0;j<
int(dominoEnergy.size());++j){
404 if (OwnerShip[j] ==
i){
405 std::vector <EcalRecHit>
temp = dominoCells[j];
406 for (
int k=0;
k<
int(temp.size());++
k){
407 dets.push_back( std::pair<DetId, float>(temp[
k].
id(), 1.) );
408 if (temp[
k].
id()==itID)
409 HasSeedCrystal =
true;
410 recHits.push_back(temp[
k]);
415 LogTrace(
"EcalClusters") <<
"Adding a cluster with: " <<
nhits;
416 LogTrace(
"EcalClusters") <<
"total E: " << LumpEnergy[
i];
417 LogTrace(
"EcalClusters") <<
"total dets: " << dets.size() ;
424 std::vector<std::pair<DetId, float> > usedHits;
425 for (
int blarg=0;blarg<
int(recHits.size());++blarg){
428 usedHits.push_back(std::make_pair<DetId, float>(recHits[blarg].
id(), 1.0));
435 if (HasSeedCrystal) {
457 if (!thisseedClusters.empty())
459 clustered_.insert(std::make_pair(clustercounter, thisseedClusters));
472 std::map<int, reco::BasicClusterCollection>::iterator mapit;
479 std::vector <reco::BasicCluster> thiscoll = mapit->second;
491 for (
int i=0;
i<
int(thiscoll.size());++
i){
493 for (
int j=0;j<
int(clustersCollection.
size());++j){
496 if (thisclus== cluster_p){
503 if (isSeed) seed = clustersCollection[j];
505 ClusterE += cluster_p.energy();
506 posX += cluster_p.energy() * cluster_p.position().X();
507 posY += cluster_p.energy() * cluster_p.position().Y();
508 posZ += cluster_p.energy() * cluster_p.position().Z();
533 SCcoll.push_back(suCl);
535 LogTrace(
"EcalClusters") <<
"Super cluster sum: " << ClusterE ;
536 LogTrace(
"EcalClusters") <<
"Made supercluster with energy E: " << suCl.
energy() ;
561 cells.push_back(SeedHit);
571 cells.push_back(UpEta);
585 cells.push_back(DownEta);
598 if (ieta2 !=
DetId(0)){
605 cells.push_back(DownEta2);
612 if (ieta1 !=
DetId(0)){
619 cells.push_back(UpEta2);
633 std::vector< std::pair<DetId, float> > dets;
638 for (
int dx = -2;
dx < 3; ++
dx)
640 for (
int dy = -2;
dy < 3; ++
dy)
646 if (thisDet !=
DetId(0))
651 dets.push_back( std::pair<DetId, float>(thisDet, 1.) );
652 energySum += hit->energy();
661 double et = energySum/cosh(pos.eta());
673 std::vector< std::pair<DetId, float> > dets;
677 for (
int dx = -2;
dx < 3; ++
dx)
679 for (
int dy = -2;
dy < 3; ++
dy)
684 if (thisDet !=
DetId(0))
689 dets.push_back( std::pair<DetId, float>(thisDet, 1.) );
697 return 1/cosh(pos.eta());
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_
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_
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
et
define resolution functions of each parameter
iterator find(key_type k)
math::XYZPoint Calculate_Location(const HitsAndFractions &iDetIds, const edm::SortedCollection< HitType > *iRecHits, const CaloSubdetectorGeometry *iSubGeom, const CaloSubdetectorGeometry *iESGeom=0)
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