#include <IslandClusterAlgo.h>

Public Types
enum	EcalPart { barrel = 0, endcap = 1 }
typedef math::XYZPoint	Point
	point in the space
enum	VerbosityLevel { pDEBUG = 0, pWARNING = 1, pINFO = 2, pERROR = 3 }
Public Member Functions
	IslandClusterAlgo ()
	IslandClusterAlgo (double ebst, double ecst, const PositionCalc &posCalc, VerbosityLevel the_verbosity=pERROR)
std::vector< reco::BasicCluster >	makeClusters (const EcalRecHitCollection hits, const CaloSubdetectorGeometry geometry, const CaloSubdetectorTopology topology_p, const CaloSubdetectorGeometry geometryES_p, EcalPart ecalPart, bool regional=false, const std::vector< EcalEtaPhiRegion > &regions=std::vector< EcalEtaPhiRegion >())
void	setVerbosity (VerbosityLevel the_verbosity)
virtual	~IslandClusterAlgo ()
Private Member Functions
void	mainSearch (const EcalRecHitCollection hits, const CaloSubdetectorGeometry geometry_p, const CaloSubdetectorTopology topology_p, const CaloSubdetectorGeometry geometryES_p, EcalPart ecalPart)
void	makeCluster (const EcalRecHitCollection hits, const CaloSubdetectorGeometry geometry_p, const CaloSubdetectorGeometry *geometryES_p)
void	searchEast (const CaloNavigator< DetId > &navigator, const CaloSubdetectorTopology *topology)
void	searchNorth (const CaloNavigator< DetId > &navigator)
void	searchSouth (const CaloNavigator< DetId > &navigator)
void	searchWest (const CaloNavigator< DetId > &navigator, const CaloSubdetectorTopology *topology)
bool	shouldBeAdded (EcalRecHitCollection::const_iterator candidate_it, EcalRecHitCollection::const_iterator previous_it)
Private Attributes
std::vector< reco::BasicCluster >	clusters_v
std::vector< std::pair< DetId, float > >	current_v
double	ecalBarrelSeedThreshold
double	ecalEndcapSeedThreshold
PositionCalc	posCalculator_
const EcalRecHitCollection *	recHits_
std::vector< EcalRecHit >	seeds
std::set< DetId >	used_s
VerbosityLevel	verbosity

Detailed Description

Definition at line 27 of file IslandClusterAlgo.h.

Member Typedef Documentation

typedef math::XYZPoint IslandClusterAlgo::Point

point in the space

Definition at line 61 of file IslandClusterAlgo.h.

Member Enumeration Documentation

enum IslandClusterAlgo::EcalPart

Enumerator:

barrel
endcap

Definition at line 31 of file IslandClusterAlgo.h.

{ barrel = 0, endcap = 1 };

enum IslandClusterAlgo::VerbosityLevel

Enumerator:

pDEBUG
pWARNING
pINFO
pERROR

Definition at line 32 of file IslandClusterAlgo.h.

{ pDEBUG = 0, pWARNING = 1, pINFO = 2, pERROR = 3 };

Constructor & Destructor Documentation

IslandClusterAlgo::IslandClusterAlgo ( ) [inline]

Definition at line 34 of file IslandClusterAlgo.h.

                      {
  }

IslandClusterAlgo::IslandClusterAlgo	(	double	ebst,
		double	ecst,
		const PositionCalc &	posCalc,
		VerbosityLevel	the_verbosity = `pERROR`
	)		`[inline]`

Definition at line 37 of file IslandClusterAlgo.h.

References posCalculator_.

                                                                                                                  : 
    ecalBarrelSeedThreshold(ebst), ecalEndcapSeedThreshold(ecst), verbosity(the_verbosity) {
    posCalculator_ = posCalc;
  }

virtual IslandClusterAlgo::~IslandClusterAlgo ( ) [inline, virtual]

Definition at line 42 of file IslandClusterAlgo.h.

    {
    }

Member Function Documentation

void IslandClusterAlgo::mainSearch	(	const EcalRecHitCollection *	hits,
		const CaloSubdetectorGeometry *	geometry_p,
		const CaloSubdetectorTopology *	topology_p,
		const CaloSubdetectorGeometry *	geometryES_p,
		EcalPart	ecalPart
	)		`[private]`

Definition at line 106 of file IslandClusterAlgo.cc.

References gather_cfg::cout, current_v, makeCluster(), pINFO, searchEast(), searchNorth(), searchSouth(), searchWest(), seeds, used_s, and verbosity.

Referenced by makeClusters().

{
  if (verbosity < pINFO)
    {
      std::cout << "Building clusters............" << std::endl;
    }

  // Loop over seeds:
  std::vector<EcalRecHit>::iterator it;
  for (it = seeds.begin(); it != seeds.end(); it++)
    {
      // make sure the current seed does not belong to a cluster already.
      if (used_s.find(it->id()) != used_s.end())
        {
          if (it == seeds.begin())
            {
              if (verbosity < pINFO)
                {
                  std::cout << "##############################################################" << std::endl;
                  std::cout << "DEBUG ALERT: Highest energy seed already belongs to a cluster!" << std::endl;
                  std::cout << "##############################################################" << std::endl;
                }
            }
          continue;
        }

      // clear the vector of hits in current cluster
      current_v.clear();

      current_v.push_back( std::pair<DetId, float>(it->id(), 1.) ); // by default hit energy fractions are set at 1.
      used_s.insert(it->id());

      // Create a navigator at the seed
      CaloNavigator<DetId> navigator(it->id(), topology_p);

      searchNorth(navigator);
      navigator.home();
      searchSouth(navigator);
      navigator.home();
      searchWest(navigator, topology_p);
      navigator.home();
      searchEast(navigator, topology_p);
 
      makeCluster(hits,geometry_p,geometryES_p);
   }
}

void IslandClusterAlgo::makeCluster	(	const EcalRecHitCollection *	hits,
		const CaloSubdetectorGeometry *	geometry_p,
		const CaloSubdetectorGeometry *	geometryES_p
	)		`[private]`

Definition at line 266 of file IslandClusterAlgo.cc.

References PositionCalc::Calculate_Location(), clusters_v, gather_cfg::cout, current_v, reco::CaloID::DET_ECAL_BARREL, reco::CaloID::DET_ECAL_ENDCAP, EcalBarrel, CaloRecHit::energy(), relval_parameters_module::energy, edm::SortedCollection< T, SORT >::find(), reco::CaloCluster::island, pINFO, posCalculator_, position, and verbosity.

Referenced by mainSearch().

{
  double energy = 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;
      if ( (*it).first.subdetId() == EcalBarrel ) {
              caloID = reco::CaloID::DET_ECAL_BARREL;
      } else {
              caloID = reco::CaloID::DET_ECAL_ENDCAP;
      }
      //      if (hit_p != 0)
      //        {
          energy += hit_p.energy();
      //        }
      //      else 
      //        {
      //          std::cout << "DEBUG ALERT: Requested rechit has gone missing from rechits map! :-S" << std::endl;
      //        }
    }

  if (verbosity < pINFO)
    { 
      std::cout << "******** NEW CLUSTER ********" << std::endl;
      std::cout << "No. of crystals = " << current_v.size() << std::endl;
      std::cout << "     Energy     = " << energy << std::endl;
      std::cout << "     Phi        = " << position.phi() << std::endl;
      std::cout << "     Eta        = " << position.eta() << std::endl;
      std::cout << "*****************************" << std::endl;
    }
  clusters_v.push_back(reco::BasicCluster(energy, position, caloID, current_v, reco::CaloCluster::island));
}

std::vector< reco::BasicCluster > IslandClusterAlgo::makeClusters	(	const EcalRecHitCollection *	hits,
		const CaloSubdetectorGeometry *	geometry,
		const CaloSubdetectorTopology *	topology_p,
		const CaloSubdetectorGeometry *	geometryES_p,
		EcalPart	ecalPart,
		bool	regional = `false`,
		const std::vector< EcalEtaPhiRegion > &	regions = `std::vector<EcalEtaPhiRegion>()`
	)

Definition at line 17 of file IslandClusterAlgo.cc.

References barrel, edm::SortedCollection< T, SORT >::begin(), clusters_v, gather_cfg::cout, ecalBarrelSeedThreshold, ecalEndcapSeedThreshold, edm::SortedCollection< T, SORT >::end(), endcap, relval_parameters_module::energy, ET, CaloSubdetectorGeometry::getGeometry(), CaloCellGeometry::getPosition(), mainSearch(), pINFO, position, recHits_, seeds, funct::sin(), python::multivaluedict::sort(), PV3DBase< T, PVType, FrameType >::theta(), dtDQMClient_cfg::threshold, used_s, and verbosity.

Referenced by EgammaHLTIslandClusterProducer::clusterizeECALPart(), IslandClusterProducer::clusterizeECALPart(), and Pi0FixedMassWindowCalibration::duringLoop().

{
  seeds.clear();
  used_s.clear();
  clusters_v.clear();

  recHits_ = hits;

  double threshold = 0;
  std::string ecalPart_string;
  if (ecalPart == endcap) 
    {
      threshold = ecalEndcapSeedThreshold;
      ecalPart_string = "EndCap";
    }
  if (ecalPart == barrel) 
    {
      threshold = ecalBarrelSeedThreshold;
      ecalPart_string = "Barrel";
    }

  if (verbosity < pINFO)
    {
      std::cout << "-------------------------------------------------------------" << std::endl;
      std::cout << "Island algorithm invoked for ECAL" << ecalPart_string << std::endl;
      std::cout << "Looking for seeds, energy threshold used = " << threshold << " GeV" <<std::endl;
    }

  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());

  if (verbosity < pINFO)
    {
      std::cout << "Total number of seeds found in event = " << seeds.size() << std::endl;
    }

  mainSearch(hits,geometry_p,topology_p,geometryES_p,ecalPart);
  sort(clusters_v.rbegin(), clusters_v.rend(), ClusterEtLess());

  if (verbosity < pINFO)
    {
      std::cout << "---------- end of main search. clusters have been sorted ----" << std::endl;
    }
  
  return clusters_v; 
}

void IslandClusterAlgo::searchEast	(	const CaloNavigator< DetId > &	navigator,
		const CaloSubdetectorTopology *	topology
	)		`[private]`

Definition at line 226 of file IslandClusterAlgo.cc.

References current_v, CaloNavigator< T >::east(), edm::SortedCollection< T, SORT >::find(), CaloNavigator< T >::home(), CaloNavigator< T >::pos(), recHits_, searchNorth(), searchSouth(), shouldBeAdded(), and used_s.

Referenced by mainSearch().

{
  DetId western = navigator.pos();
  EcalRecHitCollection::const_iterator western_it = recHits_->find(western);

  DetId eastern = navigator.east();
  if (eastern == DetId(0)) return; // This means that we went off the ECAL!
  EcalRecHitCollection::const_iterator eastern_it = recHits_->find(eastern);

  if (shouldBeAdded(eastern_it, western_it))
    {
      CaloNavigator<DetId> nsNavigator(eastern, topology);

      searchNorth(nsNavigator);
      nsNavigator.home();
      searchSouth(nsNavigator);
      nsNavigator.home();
      searchEast(navigator, topology);

      current_v.push_back( std::pair<DetId, float>(eastern, 1.)); // by default hit energy fractions are set at 1.
      used_s.insert(eastern);
    }
}

void IslandClusterAlgo::searchNorth ( const CaloNavigator< DetId > & navigator ) [private]

Definition at line 158 of file IslandClusterAlgo.cc.

References current_v, edm::SortedCollection< T, SORT >::find(), CaloNavigator< T >::north(), CaloNavigator< T >::pos(), recHits_, shouldBeAdded(), and used_s.

Referenced by mainSearch(), searchEast(), and searchWest().

{
  DetId southern = navigator.pos();

  DetId northern = navigator.north();
  if (northern == DetId(0)) return; // This means that we went off the ECAL!
  // if the crystal to the north belongs to another cluster return
  if (used_s.find(northern) != used_s.end()) return;


  EcalRecHitCollection::const_iterator southern_it = recHits_->find(southern);
  EcalRecHitCollection::const_iterator northern_it = recHits_->find(northern);

  if (shouldBeAdded(northern_it, southern_it))
    {
      current_v.push_back( std::pair<DetId, float>(northern, 1.)); // by default hit energy fractions are set at 1.
      used_s.insert(northern);
      searchNorth(navigator);
    }
}

void IslandClusterAlgo::searchSouth ( const CaloNavigator< DetId > & navigator ) [private]

Definition at line 180 of file IslandClusterAlgo.cc.

References current_v, edm::SortedCollection< T, SORT >::find(), CaloNavigator< T >::pos(), recHits_, shouldBeAdded(), CaloNavigator< T >::south(), and used_s.

Referenced by mainSearch(), searchEast(), and searchWest().

{
  DetId northern = navigator.pos();

  DetId southern = navigator.south();
  if (southern == DetId(0)) return; // This means that we went off the ECAL!
  if (used_s.find(southern) != used_s.end()) return;


  EcalRecHitCollection::const_iterator northern_it = recHits_->find(northern);
  EcalRecHitCollection::const_iterator southern_it = recHits_->find(southern);

  if (shouldBeAdded(southern_it, northern_it))
    {
      current_v.push_back( std::pair<DetId, float>(southern, 1.)); // by default hit energy fractions are set at 1.
      used_s.insert(southern);
      searchSouth(navigator);
    }
}

void IslandClusterAlgo::searchWest	(	const CaloNavigator< DetId > &	navigator,
		const CaloSubdetectorTopology *	topology
	)		`[private]`

Definition at line 201 of file IslandClusterAlgo.cc.

References current_v, edm::SortedCollection< T, SORT >::find(), CaloNavigator< T >::home(), CaloNavigator< T >::pos(), recHits_, searchNorth(), searchSouth(), shouldBeAdded(), used_s, and CaloNavigator< T >::west().

Referenced by mainSearch().

{
  DetId eastern = navigator.pos();
  EcalRecHitCollection::const_iterator eastern_it = recHits_->find(eastern);

  DetId western = navigator.west();
  if (western == DetId(0)) return; // This means that we went off the ECAL!
  EcalRecHitCollection::const_iterator western_it = recHits_->find(western);

  if (shouldBeAdded(western_it, eastern_it))
    {
      CaloNavigator<DetId> nsNavigator(western, topology);

      searchNorth(nsNavigator);
      nsNavigator.home();
      searchSouth(nsNavigator);
      nsNavigator.home();
      searchWest(navigator, topology);

      current_v.push_back( std::pair<DetId, float>(western, 1.)); // by default hit energy fractions are set at 1.
      used_s.insert(western);
    }
}

void IslandClusterAlgo::setVerbosity ( VerbosityLevel the_verbosity ) [inline]

Definition at line 46 of file IslandClusterAlgo.h.

References verbosity.

    {
      verbosity = the_verbosity;
    }

bool IslandClusterAlgo::shouldBeAdded	(	EcalRecHitCollection::const_iterator	candidate_it,
		EcalRecHitCollection::const_iterator	previous_it
	)		`[private]`

Definition at line 252 of file IslandClusterAlgo.cc.

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

Referenced by searchEast(), searchNorth(), searchSouth(), and searchWest().

{
  // crystal should not be included...
  if ( (candidate_it == recHits_->end())                 || // ...if it does not correspond to a hit
       (used_s.find(candidate_it->id()) != used_s.end()) || // ...if it already belongs to a cluster
       (candidate_it->energy() <= 0)                     || // ...if it has a negative or zero energy
       (candidate_it->energy() > previous_it->energy()))    // ...or if the previous crystal had lower E
    {
      return false;
    }
  return true;
}