#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< EcalEtaPhiRegion > &regions=std::vector< EcalEtaPhiRegion >())

	Multi5x5ClusterAlgo ()

	Multi5x5ClusterAlgo (double ebst, double ecst, 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 28 of file Multi5x5ClusterAlgo.h.

Member Typedef Documentation

typedef math::XYZPoint Multi5x5ClusterAlgo::Point

point in the space

Definition at line 81 of file Multi5x5ClusterAlgo.h.

Constructor & Destructor Documentation

Multi5x5ClusterAlgo::Multi5x5ClusterAlgo ( )

inline

Definition at line 57 of file Multi5x5ClusterAlgo.h.

57 {

58 }

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

inline

Definition at line 60 of file Multi5x5ClusterAlgo.h.

References posCalculator_, python.multivaluedict::sort(), 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 66 of file Multi5x5ClusterAlgo.h.

67 {

68 }

Member Function Documentation

void Multi5x5ClusterAlgo::addCrystal ( const DetId & det )

private

Definition at line 407 of file Multi5x5ClusterAlgo.cc.

References current_v, edm::SortedCollection< T, SORT >::end(), edm::SortedCollection< T, SORT >::find(), recHits_, and used_s.

Referenced by prepareCluster().

 {   
 
     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 298 of file Multi5x5ClusterAlgo.cc.

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

Referenced by mainSearch().

 {
 
     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();
 
     std::vector<DetId>::const_iterator detItr;
     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 136 of file Multi5x5ClusterAlgo.cc.

References PositionCalc::Calculate_Location(), canSeed_s, checkMaxima(), clusters_v, current_v, detector_, Multi5x5ClusterAlgo::ProtoBasicCluster::energy(), spr::find(), edm::SortedCollection< T, SORT >::find(), Multi5x5ClusterAlgo::ProtoBasicCluster::hits(), EcalRecHit::id(), LogTrace, makeCluster(), reco::CaloCluster::multi5x5, posCalculator_, position, prepareCluster(), protoClusters_, reassignSeedCrysToClusterItSeeds_, query::result, Multi5x5ClusterAlgo::ProtoBasicCluster::seed(), seeds, python.multivaluedict::sort(), used_s, v_chstatus_, and whichClusCrysBelongsTo_.

Referenced by makeClusters().

 {
 
     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 (recoFlag based)
         uint32_t rhFlag = (*it).recoFlag();
         std::vector<int>::const_iterator vit = std::find( v_chstatus_.begin(), v_chstatus_.end(), rhFlag );
         if ( vit != v_chstatus_.end() ) continue; // the recHit has to be excluded from seeding
 
         // 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.size() > 0) 
         {
             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 237 of file Multi5x5ClusterAlgo.cc.

References PositionCalc::Calculate_Location(), current_v, reco::CaloID::DET_ECAL_BARREL, reco::CaloID::DET_ECAL_ENDCAP, EcalBarrel, CaloRecHit::energy(), relval_parameters_module::energy, edm::SortedCollection< T, SORT >::find(), LogTrace, posCalculator_, position, protoClusters_, reassignSeedCrysToClusterItSeeds_, used_s, and whichClusCrysBelongsTo_.

Referenced by mainSearch().

 {
 
     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< EcalEtaPhiRegion > &	regions = `std::vector<EcalEtaPhiRegion>()`
	)

Definition at line 46 of file Multi5x5ClusterAlgo.cc.

References edm::SortedCollection< T, SORT >::begin(), canSeed_s, clusters_v, reco::CaloID::DET_ECAL_BARREL, reco::CaloID::DET_ECAL_ENDCAP, reco::CaloID::DET_NONE, detector_, ecalBarrelSeedThreshold, ecalEndcapSeedThreshold, edm::SortedCollection< T, SORT >::end(), relval_parameters_module::energy, ET, CaloSubdetectorGeometry::getGeometry(), CaloCellGeometry::getPosition(), LogTrace, mainSearch(), position, recHits_, seeds, funct::sin(), python.multivaluedict::sort(), PV3DBase< T, PVType, FrameType >::theta(), dtDQMClient_cfg::threshold, and used_s.

Referenced by Multi5x5ClusterProducer::clusterizeECALPart(), and EgammaHLTMulti5x5ClusterProducer::clusterizeECALPart().

 {
     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!
 
             const CaloCellGeometry *thisCell = geometry_p->getGeometry(it->id());
             GlobalPoint position = thisCell->getPosition();
 
             // Require that RecHit is within clustering region in case
             // of regional reconstruction
             bool withinRegion = false;
             if (regional) {
                 std::vector<EcalEtaPhiRegion>::const_iterator region;
                 for (region=regions.begin(); region!=regions.end(); region++) {
                     if (region->inRegion(position)) {
                         withinRegion =  true;
                         break;
                     }
                 }
             }
 
             if (!regional || withinRegion) {
                 float ET = it->energy() * sin(position.theta());
                 if (ET > threshold) seeds.push_back(*it);
             }
         }
 
     }
 
     sort(seeds.begin(), seeds.end(), EcalRecHitLess());
 
     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(), ClusterEtLess());
 
     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 348 of file Multi5x5ClusterAlgo.cc.

References abs, addCrystal(), canSeed_s, CaloNavigator< T >::home(), and CaloNavigator< T >::offsetBy().

Referenced by mainSearch().

 {
 
     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;
 }