#include <Multi5x5ClusterAlgo.h>

Classes
class	ProtoBasicCluster
Public Types
typedef math::XYZPoint	Point
	point in the space
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 (double ebst, double ecst, std::vector< int > v_chstatus, const PositionCalc &posCalc, bool reassignSeedCrysToClusterItSeeds=false)
	Multi5x5ClusterAlgo ()
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.
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.

                        {
        }

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 ( ) [inline, virtual]

Definition at line 66 of file Multi5x5ClusterAlgo.h.

        {
        }

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 edm::SortedCollection< T, SORT >::end(), edm::SortedCollection< T, SORT >::find(), spr::find(), i, recHits_, and v_chstatus_.

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(), edm::SortedCollection< T, SORT >::find(), Multi5x5ClusterAlgo::ProtoBasicCluster::hits(), EcalRecHit::id(), EcalRecHit::kGood, 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         
        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.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(), AlCaHLTBitMon_QueryRunRegistry::string, PV3DBase< T, PVType, FrameType >::theta(), dtDQMClient_cfg::threshold, and used_s.

Referenced by EgammaHLTMulti5x5ClusterProducer::clusterizeECALPart(), and Multi5x5ClusterProducer::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(), and canSeed_s.

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