#include <Multi5x5ClusterAlgo.h>

Classes
class	ProtoBasicCluster

Public Types
typedef math::XYZPoint	Point
	point in the space More...

Public Member Functions
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 >())

	Multi5x5ClusterAlgo ()

	Multi5x5ClusterAlgo (double ebst, double ecst, const std::vector< int > &v_chstatus, const PositionCalc &posCalc, bool reassignSeedCrysToClusterItSeeds=false)

virtual	~Multi5x5ClusterAlgo ()

Private Member Functions
void	addCrystal (const DetId &det)

bool	checkMaxima (CaloNavigator< DetId > &navigator, const EcalRecHitCollection *hits)

void	mainSearch (const EcalRecHitCollection hits, const CaloSubdetectorGeometry geometry_p, const CaloSubdetectorTopology topology_p, const CaloSubdetectorGeometry geometryES_p)

void	makeCluster (const EcalRecHitCollection hits, const CaloSubdetectorGeometry geometry_p, const CaloSubdetectorGeometry *geometryES_p, const EcalRecHitCollection::const_iterator &seedIt, bool seedOutside)

void	prepareCluster (CaloNavigator< DetId > &navigator, const EcalRecHitCollection hits, const CaloSubdetectorGeometry geometry)

Private Attributes
std::set< DetId >	canSeed_s

std::vector< reco::BasicCluster >	clusters_v

std::vector< std::pair< DetId, float > >	current_v

reco::CaloID::Detectors	detector_
	The ecal region used. More...

double	ecalBarrelSeedThreshold

double	ecalEndcapSeedThreshold

PositionCalc	posCalculator_

std::vector< ProtoBasicCluster >	protoClusters_

bool	reassignSeedCrysToClusterItSeeds_

const EcalRecHitCollection *	recHits_

std::vector< EcalRecHit >	seeds

std::set< DetId >	used_s

std::vector< int >	v_chstatus_

std::vector< std::pair< DetId, int > >	whichClusCrysBelongsTo_

Detailed Description

Definition at line 27 of file Multi5x5ClusterAlgo.h.

Member Typedef Documentation

typedef math::XYZPoint Multi5x5ClusterAlgo::Point

point in the space

Definition at line 87 of file Multi5x5ClusterAlgo.h.

Constructor & Destructor Documentation

Multi5x5ClusterAlgo::Multi5x5ClusterAlgo ( )

inline

Definition at line 59 of file Multi5x5ClusterAlgo.h.

59 {}

Multi5x5ClusterAlgo::Multi5x5ClusterAlgo	(	double	ebst,
		double	ecst,
		const std::vector< int > &	v_chstatus,
		const PositionCalc &	posCalc,
		bool	reassignSeedCrysToClusterItSeeds = `false`
	)

inline

Definition at line 61 of file Multi5x5ClusterAlgo.h.

References posCalculator_, and v_chstatus_.

       : ecalBarrelSeedThreshold(ebst),
         ecalEndcapSeedThreshold(ecst),
         v_chstatus_(v_chstatus),
         reassignSeedCrysToClusterItSeeds_(reassignSeedCrysToClusterItSeeds) {
     posCalculator_ = posCalc;
     std::sort(v_chstatus_.begin(), v_chstatus_.end());
   }

virtual Multi5x5ClusterAlgo::~Multi5x5ClusterAlgo ( )

inlinevirtual

Definition at line 74 of file Multi5x5ClusterAlgo.h.

References hgcalTestNeighbor_cfi::detector, Multi5x5ClusterAlgo::ProtoBasicCluster::hits(), and makeClusters().

74 {}

Member Function Documentation

void Multi5x5ClusterAlgo::addCrystal ( const DetId & det )

private

Definition at line 376 of file Multi5x5ClusterAlgo.cc.

                                                      {
   EcalRecHitCollection::const_iterator thisIt = recHits_->find(det);
   if ((thisIt != recHits_->end()) && (thisIt->id() != DetId(0))) {
     if ((used_s.find(thisIt->id()) == used_s.end())) {
       //std::cout << "   ... this is a good crystal and will be added" << std::endl;
       current_v.push_back(std::pair<DetId, float>(det, 1.));  // by default hit energy fractions are set at 1.
       used_s.insert(det);
     }
   }
 }

bool Multi5x5ClusterAlgo::checkMaxima	(	CaloNavigator< DetId > &	navigator,
		const EcalRecHitCollection *	hits
	)

private

Definition at line 283 of file Multi5x5ClusterAlgo.cc.

References CaloNavigator< T, TOPO >::east(), spr::find(), edm::SortedCollection< T, SORT >::find(), CaloNavigator< T, TOPO >::home(), mps_fire::i, CaloNavigator< T, TOPO >::north(), CaloNavigator< T, TOPO >::pos(), CaloNavigator< T, TOPO >::south(), and CaloNavigator< T, TOPO >::west().

                                                                                                        {
   bool maxima = true;
   EcalRecHitCollection::const_iterator thisHit;
   EcalRecHitCollection::const_iterator seedHit = hits->find(navigator.pos());
   double seedEnergy = seedHit->energy();
 
   std::vector<DetId> swissCrossVec;
   swissCrossVec.clear();
 
   swissCrossVec.push_back(navigator.west());
   navigator.home();
   swissCrossVec.push_back(navigator.east());
   navigator.home();
   swissCrossVec.push_back(navigator.north());
   navigator.home();
   swissCrossVec.push_back(navigator.south());
   navigator.home();
 
   for (unsigned int i = 0; i < swissCrossVec.size(); ++i) {
     // look for this hit
     thisHit = recHits_->find(swissCrossVec[i]);
 
     // continue if this hit was not found
     if ((swissCrossVec[i] == DetId(0)) || thisHit == recHits_->end())
       continue;
 
     // the recHit has to be skipped in the local maximum search if it was found
     // in the map of channels to be excluded
     uint32_t rhFlag = thisHit->recoFlag();
     std::vector<int>::const_iterator vit = std::find(v_chstatus_.begin(), v_chstatus_.end(), rhFlag);
     if (vit != v_chstatus_.end())
       continue;
 
     // if this crystal has more energy than the seed then we do
     // not have a local maxima
     if (thisHit->energy() > seedEnergy) {
       maxima = false;
       break;
     }
   }
 
   return maxima;
 }

void Multi5x5ClusterAlgo::mainSearch	(	const EcalRecHitCollection *	hits,
		const CaloSubdetectorGeometry *	geometry_p,
		const CaloSubdetectorTopology *	topology_p,
		const CaloSubdetectorGeometry *	geometryES_p
	)

private

Definition at line 133 of file Multi5x5ClusterAlgo.cc.

References Multi5x5ClusterAlgo::ProtoBasicCluster::energy(), edm::SortedCollection< T, SORT >::find(), hfClusterShapes_cfi::hits, Multi5x5ClusterAlgo::ProtoBasicCluster::hits(), EcalRecHit::id(), EcalRecHit::kGood, LogTrace, reco::CaloCluster::multi5x5, HLT_2018_cff::navigator, position, mps_fire::result, and Multi5x5ClusterAlgo::ProtoBasicCluster::seed().

                                                                                   {
   LogTrace("EcalClusters") << "Building clusters............";
 
   // Loop over seeds:
   std::vector<EcalRecHit>::iterator it;
   for (it = seeds.begin(); it != seeds.end(); it++) {
     // check if this crystal is able to seed
     // (event though it is already used)
     bool usedButCanSeed = false;
     if (canSeed_s.find(it->id()) != canSeed_s.end())
       usedButCanSeed = true;
 
     // avoid seeding for anomalous channels
     if (!it->checkFlag(EcalRecHit::kGood)) {  // if rechit is good, no need for further checks
       if (it->checkFlags(v_chstatus_))
         continue;
     }
 
     // make sure the current seed does not belong to a cluster already.
     if ((used_s.find(it->id()) != used_s.end()) && (usedButCanSeed == false)) {
       if (it == seeds.begin()) {
         LogTrace("EcalClusters") << "##############################################################";
         LogTrace("EcalClusters") << "DEBUG ALERT: Highest energy seed already belongs to a cluster!";
         LogTrace("EcalClusters") << "##############################################################";
       }
 
       // seed crystal is used or is used and cannot seed a cluster
       // so continue to the next seed crystal...
       continue;
     }
 
     // clear the vector of hits in current cluster
     current_v.clear();
 
     // Create a navigator at the seed and get seed
     // energy
     CaloNavigator<DetId> navigator(it->id(), topology_p);
     DetId seedId = navigator.pos();
     EcalRecHitCollection::const_iterator seedIt = hits->find(seedId);
     navigator.setHome(seedId);
 
     // Is the seed a local maximum?
     bool localMaxima = checkMaxima(navigator, hits);
 
     if (localMaxima) {
       // build the 5x5 taking care over which crystals
       // can seed new clusters and which can't
       prepareCluster(navigator, hits, geometry_p);
     }
 
     // If some crystals in the current vector then
     // make them into a cluster
     if (!current_v.empty()) {
       makeCluster(hits, geometry_p, geometryES_p, seedIt, usedButCanSeed);
     }
 
   }  // End loop on seed crystals
 
   if (reassignSeedCrysToClusterItSeeds_) {
     std::sort(whichClusCrysBelongsTo_.begin(), whichClusCrysBelongsTo_.end(), PairSortByFirst<DetId, int>());
 
     for (size_t clusNr = 0; clusNr < protoClusters_.size(); clusNr++) {
       if (!protoClusters_[clusNr].containsSeed()) {
         const EcalRecHit& seedHit = protoClusters_[clusNr].seed();
         typedef std::vector<std::pair<DetId, int> >::iterator It;
         std::pair<It, It> result = std::equal_range(whichClusCrysBelongsTo_.begin(),
                                                     whichClusCrysBelongsTo_.end(),
                                                     seedHit.id(),
                                                     PairSortByFirst<DetId, int>());
 
         if (result.first != result.second)
           protoClusters_[result.first->second].removeHit(seedHit);
         protoClusters_[clusNr].addSeed();
       }
     }
   }
 
   for (size_t clusNr = 0; clusNr < protoClusters_.size(); clusNr++) {
     const ProtoBasicCluster& protoCluster = protoClusters_[clusNr];
     Point position;
     position = posCalculator_.Calculate_Location(protoCluster.hits(), hits, geometry_p, geometryES_p);
     clusters_v.push_back(reco::BasicCluster(protoCluster.energy(),
                                             position,
                                             reco::CaloID(detector_),
                                             protoCluster.hits(),
                                             reco::CaloCluster::multi5x5,
                                             protoCluster.seed().id()));
   }
 
   protoClusters_.clear();
   whichClusCrysBelongsTo_.clear();
 }

void Multi5x5ClusterAlgo::makeCluster	(	const EcalRecHitCollection *	hits,
		const CaloSubdetectorGeometry *	geometry_p,
		const CaloSubdetectorGeometry *	geometryES_p,
		const EcalRecHitCollection::const_iterator &	seedIt,
		bool	seedOutside
	)

private

Definition at line 229 of file Multi5x5ClusterAlgo.cc.

References reco::CaloID::DET_ECAL_BARREL, reco::CaloID::DET_ECAL_ENDCAP, EcalBarrel, HCALHighEnergyHPDFilter_cfi::energy, EcalRecHit::energy(), edm::SortedCollection< T, SORT >::find(), LogTrace, and position.

                                                         {
   double energy = 0;
   //double chi2   = 0;
   reco::CaloID caloID;
   Point position;
   position = posCalculator_.Calculate_Location(current_v, hits, geometry, geometryES);
 
   std::vector<std::pair<DetId, float> >::iterator it;
   for (it = current_v.begin(); it != current_v.end(); it++) {
     EcalRecHitCollection::const_iterator itt = hits->find((*it).first);
     EcalRecHit hit_p = *itt;
     energy += hit_p.energy();
     //chi2 += 0;
     if ((*it).first.subdetId() == EcalBarrel) {
       caloID = reco::CaloID::DET_ECAL_BARREL;
     } else {
       caloID = reco::CaloID::DET_ECAL_ENDCAP;
     }
   }
   //chi2 /= energy;
 
   LogTrace("EcalClusters") << "******** NEW CLUSTER ********";
   LogTrace("EcalClusters") << "No. of crystals = " << current_v.size();
   LogTrace("EcalClusters") << "     Energy     = " << energy;
   LogTrace("EcalClusters") << "     Phi        = " << position.phi();
   LogTrace("EcalClusters") << "     Eta " << position.eta();
   LogTrace("EcalClusters") << "*****************************";
 
   // to be a valid cluster the cluster energy
   // must be at least the seed energy
   double seedEnergy = seedIt->energy();
   if ((seedOutside && energy >= 0) || (!seedOutside && energy >= seedEnergy)) {
     if (reassignSeedCrysToClusterItSeeds_) {  //if we're not doing this, we dont need this info so lets not bother filling it
       for (size_t hitNr = 0; hitNr < current_v.size(); hitNr++)
         whichClusCrysBelongsTo_.push_back(std::pair<DetId, int>(current_v[hitNr].first, protoClusters_.size()));
     }
     protoClusters_.push_back(ProtoBasicCluster(energy, *seedIt, current_v));
 
     // clusters_v.push_back(reco::BasicCluster(energy, position, reco::CaloID(detector_), current_v, reco::CaloCluster::multi5x5, seedIt->id()));
 
     // if no valid cluster was built,
     // then free up these crystals to be used in the next...
   } else {
     std::vector<std::pair<DetId, float> >::iterator iter;
     for (iter = current_v.begin(); iter != current_v.end(); iter++) {
       used_s.erase(iter->first);
     }  //for(iter)
   }    //else
 }

std::vector< reco::BasicCluster > Multi5x5ClusterAlgo::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>()`
	)

Definition at line 54 of file Multi5x5ClusterAlgo.cc.

References edm::SortedCollection< T, SORT >::begin(), reco::CaloID::DET_ECAL_BARREL, reco::CaloID::DET_ECAL_ENDCAP, reco::CaloID::DET_NONE, edm::SortedCollection< T, SORT >::end(), HCALHighEnergyHPDFilter_cfi::energy, ET, CaloSubdetectorGeometry::getGeometry(), hfClusterShapes_cfi::hits, isClusterEtLess(), LogTrace, HLT_2018_cff::region, AlCaHLTBitMon_QueryRunRegistry::string, MessageLogger_cff::threshold, x, and y.

Referenced by Multi5x5ClusterProducer::clusterizeECALPart(), and ~Multi5x5ClusterAlgo().

                                                                                                                      {
   seeds.clear();
   used_s.clear();
   canSeed_s.clear();
   clusters_v.clear();
 
   recHits_ = hits;
 
   double threshold = 0;
   std::string ecalPart_string;
   detector_ = reco::CaloID::DET_NONE;
   if (detector == reco::CaloID::DET_ECAL_ENDCAP) {
     detector_ = reco::CaloID::DET_ECAL_ENDCAP;
     threshold = ecalEndcapSeedThreshold;
     ecalPart_string = "EndCap";
   }
   if (detector == reco::CaloID::DET_ECAL_BARREL) {
     detector_ = reco::CaloID::DET_ECAL_BARREL;
     threshold = ecalBarrelSeedThreshold;
     ecalPart_string = "Barrel";
   }
 
   LogTrace("EcalClusters") << "-------------------------------------------------------------";
   LogTrace("EcalClusters") << "Island algorithm invoked for ECAL" << ecalPart_string;
   LogTrace("EcalClusters") << "Looking for seeds, energy threshold used = " << threshold << " GeV";
 
   int nregions = 0;
   if (regional)
     nregions = regions.size();
 
   if (!regional || nregions) {
     EcalRecHitCollection::const_iterator it;
     for (it = hits->begin(); it != hits->end(); it++) {
       double energy = it->energy();
       if (energy < threshold)
         continue;  // need to check to see if this line is useful!
 
       auto thisCell = geometry_p->getGeometry(it->id());
       // Require that RecHit is within clustering region in case
       // of regional reconstruction
       bool withinRegion = false;
       if (regional) {
         std::vector<RectangularEtaPhiRegion>::const_iterator region;
         for (region = regions.begin(); region != regions.end(); region++) {
           if (region->inRegion(thisCell->etaPos(), thisCell->phiPos())) {
             withinRegion = true;
             break;
           }
         }
       }
 
       if (!regional || withinRegion) {
         float ET = it->energy() * thisCell->getPosition().basicVector().unit().perp();
         if (ET > threshold)
           seeds.push_back(*it);
       }
     }
   }
 
   sort(seeds.begin(), seeds.end(), [](auto const& x, auto const& y) { return (x.energy() > y.energy()); });
 
   LogTrace("EcalClusters") << "Total number of seeds found in event = " << seeds.size();
 
   mainSearch(hits, geometry_p, topology_p, geometryES_p);
   sort(clusters_v.rbegin(), clusters_v.rend(), isClusterEtLess);
 
   LogTrace("EcalClusters") << "---------- end of main search. clusters have been sorted ----";
 
   return clusters_v;
 }

void Multi5x5ClusterAlgo::prepareCluster	(	CaloNavigator< DetId > &	navigator,
		const EcalRecHitCollection *	hits,
		const CaloSubdetectorGeometry *	geometry
	)

private

Definition at line 327 of file Multi5x5ClusterAlgo.cc.

References funct::abs(), PVValHelper::dx, PVValHelper::dy, CaloNavigator< T, TOPO >::home(), and CaloNavigator< T, TOPO >::offsetBy().

                                                                                   {
   DetId thisDet;
   std::set<DetId>::iterator setItr;
 
   // now add the 5x5 taking care to mark the edges
   // as able to seed and where overlapping in the central
   // region with crystals that were previously able to seed
   // change their status so they are not able to seed
   //std::cout << std::endl;
   for (int dx = -2; dx < 3; ++dx) {
     for (int dy = -2; dy < 3; ++dy) {
       // navigate in free steps forming
       // a full 5x5
       thisDet = navigator.offsetBy(dx, dy);
       navigator.home();
 
       // add the current crystal
       //std::cout << "adding " << dx << ", " << dy << std::endl;
       addCrystal(thisDet);
 
       // now consider if we are in an edge (outer 16)
       // or central (inner 9) region
       if ((abs(dx) > 1) || (abs(dy) > 1)) {
         // this is an "edge" so should be allowed to seed
         // provided it is not already used
         //std::cout << "   setting can seed" << std::endl;
         canSeed_s.insert(thisDet);
       }  // end if "edge"
       else {
         // or else we are in the central 3x3
         // and must remove any of these crystals from the canSeed set
         setItr = canSeed_s.find(thisDet);
         if (setItr != canSeed_s.end()) {
           //std::cout << "   unsetting can seed" << std::endl;
           canSeed_s.erase(setItr);
         }
       }  // end if "centre"
 
     }  // end loop on dy
 
   }  // end loop on dx
 
   //std::cout << "*** " << std::endl;
   //std::cout << " current_v contains " << current_v.size() << std::endl;
   //std::cout << "*** " << std::endl;
 }

Member Data Documentation

std::set<DetId> Multi5x5ClusterAlgo::canSeed_s

private

Definition at line 110 of file Multi5x5ClusterAlgo.h.

std::vector<reco::BasicCluster> Multi5x5ClusterAlgo::clusters_v

private

Definition at line 117 of file Multi5x5ClusterAlgo.h.

std::vector<std::pair<DetId, float> > Multi5x5ClusterAlgo::current_v

private

Definition at line 114 of file Multi5x5ClusterAlgo.h.

reco::CaloID::Detectors Multi5x5ClusterAlgo::detector_

private

The ecal region used.

Definition at line 94 of file Multi5x5ClusterAlgo.h.

double Multi5x5ClusterAlgo::ecalBarrelSeedThreshold

private

Definition at line 97 of file Multi5x5ClusterAlgo.h.

double Multi5x5ClusterAlgo::ecalEndcapSeedThreshold

private

Definition at line 98 of file Multi5x5ClusterAlgo.h.

PositionCalc Multi5x5ClusterAlgo::posCalculator_

private

Definition at line 91 of file Multi5x5ClusterAlgo.h.

Referenced by Multi5x5ClusterAlgo().

std::vector<ProtoBasicCluster> Multi5x5ClusterAlgo::protoClusters_

private

Definition at line 118 of file Multi5x5ClusterAlgo.h.

bool Multi5x5ClusterAlgo::reassignSeedCrysToClusterItSeeds_

private

Definition at line 122 of file Multi5x5ClusterAlgo.h.

const EcalRecHitCollection* Multi5x5ClusterAlgo::recHits_

private

Definition at line 101 of file Multi5x5ClusterAlgo.h.

std::vector<EcalRecHit> Multi5x5ClusterAlgo::seeds

private

Definition at line 104 of file Multi5x5ClusterAlgo.h.

std::set<DetId> Multi5x5ClusterAlgo::used_s

private

Definition at line 109 of file Multi5x5ClusterAlgo.h.

std::vector<int> Multi5x5ClusterAlgo::v_chstatus_

private

Definition at line 120 of file Multi5x5ClusterAlgo.h.

Referenced by Multi5x5ClusterAlgo().

std::vector<std::pair<DetId, int> > Multi5x5ClusterAlgo::whichClusCrysBelongsTo_

private

Definition at line 106 of file Multi5x5ClusterAlgo.h.

Classes

Public Types

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Member Typedef Documentation

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation