#include <HybridClusterAlgo.h>

Public Member Functions
	HybridClusterAlgo ()

	HybridClusterAlgo (double eb_str, int step, double ethres, double eseed, double xi, bool useEtForXi, double ewing, const std::vector< int > &v_chstatus, const PositionCalc &posCalculator, bool dynamicEThres=false, double eThresA=0, double eThresB=0.1, const std::vector< int > &severityToExclude=std::vector< int >().operator=(std::vector< int >(999)), bool excludeFromCluster=false)

void	mainSearch (const EcalRecHitCollection hits, const CaloSubdetectorGeometry geometry)

void	makeClusters (const EcalRecHitCollection , const CaloSubdetectorGeometry geometry, reco::BasicClusterCollection &basicClusters, const EcalSeverityLevelAlgo *sevLv, bool regional=false, const std::vector< EcalEtaPhiRegion > &regions=std::vector< EcalEtaPhiRegion >())

double	makeDomino (EcalBarrelNavigatorHT &navigator, std::vector< EcalRecHit > &cells)

reco::SuperClusterCollection	makeSuperClusters (const reco::CaloClusterPtrVector &)

void	setDynamicPhiRoad (const edm::ParameterSet &bremRecoveryPset)

	~HybridClusterAlgo ()

Private Types
typedef math::XYZPoint	Point

Private Member Functions
double	e2Et (EcalBarrelNavigatorHT &navigator, const EcalRecHitCollection hits, const CaloSubdetectorGeometry geometry)

double	et25 (EcalBarrelNavigatorHT &navigator, const EcalRecHitCollection hits, const CaloSubdetectorGeometry geometry)

Private Attributes
std::map< int, std::vector < reco::BasicCluster > >	clustered_

bool	dynamicEThres_

bool	dynamicPhiRoad_

double	eb_st

double	Eseed

double	eThres_

double	eThresA_

double	eThresB_

double	Ewing

std::set< DetId >	excludedCrys_

bool	excludeFromCluster_

BremRecoveryPhiRoadAlgo *	phiRoadAlgo_

int	phiSteps_

PositionCalc	posCalculator_

const EcalRecHitCollection *	recHits_

std::vector< reco::BasicCluster >	seedClus_

std::vector< EcalRecHit >	seeds

float	severityRecHitThreshold_

float	severitySpikeThreshold_

EcalBarrelHardcodedTopology *	topo_

std::set< DetId >	useddetids

bool	UseEtForXi

std::vector< int >	v_chstatus_

std::vector< int >	v_severitylevel_

double	Xi

Detailed Description

Definition at line 33 of file HybridClusterAlgo.h.

Member Typedef Documentation

typedef math::XYZPoint HybridClusterAlgo::Point

private

Definition at line 37 of file HybridClusterAlgo.h.

Constructor & Destructor Documentation

HybridClusterAlgo::HybridClusterAlgo ( )

inline

Definition at line 123 of file HybridClusterAlgo.h.

123 { }

HybridClusterAlgo::HybridClusterAlgo	(	double	eb_str,
		int	step,
		double	ethres,
		double	eseed,
		double	xi,
		bool	useEtForXi,
		double	ewing,
		const std::vector< int > &	v_chstatus,
		const PositionCalc &	posCalculator,
		bool	dynamicEThres = `false`,
		double	eThresA = `0`,
		double	eThresB = `0.1`,
		const std::vector< int > &	severityToExclude = `std::vector<int>().operator=(std::vector<int>(999))`,
		bool	excludeFromCluster = `false`
	)

Definition at line 15 of file HybridClusterAlgo.cc.

References dynamicEThres_, dynamicPhiRoad_, Eseed, eThresA_, eThresB_, LogTrace, posCalculator_, python.multivaluedict::sort(), topo_, UseEtForXi, v_chstatus_, and Xi.

                        :
   
   eb_st(eb_str), phiSteps_(step), 
   eThres_(ethres), eThresA_(eThresA), eThresB_(eThresB),
   Eseed(eseed), Xi(xi), UseEtForXi(useEtForXi), Ewing(ewing), 
   dynamicEThres_(dynamicEThres),
   v_chstatus_(v_chstatus), v_severitylevel_(severityToExclude),
   //severityRecHitThreshold_(severityRecHitThreshold),
   //severitySpikeThreshold_(severitySpikeThreshold),
   excludeFromCluster_(excludeFromCluster)
   
   
 {
 
   dynamicPhiRoad_ = false;
   LogTrace("EcalClusters") << "dynamicEThres: " << dynamicEThres_ << " : A,B " << eThresA_ << ", " << eThresB_;
   LogTrace("EcalClusters") << "Eseed: " << Eseed << " , Xi: " << Xi << ", UseEtForXi: " << UseEtForXi;
   
   //if (dynamicPhiRoad_) phiRoadAlgo_ = new BremRecoveryPhiRoadAlgo(bremRecoveryPset);
   posCalculator_ = posCalculator;
   topo_ = new EcalBarrelHardcodedTopology();
   
   std::sort( v_chstatus_.begin(), v_chstatus_.end() );
   
   
 }

HybridClusterAlgo::~HybridClusterAlgo ( )

inline

Definition at line 147 of file HybridClusterAlgo.h.

References dynamicPhiRoad_, phiRoadAlgo_, and topo_.

   {
      if (dynamicPhiRoad_) delete phiRoadAlgo_;
      delete topo_; 
     //     delete SCShape_;
   } 

Member Function Documentation

double HybridClusterAlgo::e2Et	(	EcalBarrelNavigatorHT &	navigator,
		const EcalRecHitCollection *	hits,
		const CaloSubdetectorGeometry *	geometry
	)

private

Definition at line 666 of file HybridClusterAlgo.cc.

References PositionCalc::Calculate_Location(), edm::SortedCollection< T, SORT >::end(), edm::SortedCollection< T, SORT >::find(), CaloNavigator< T, TOPO >::home(), CaloNavigator< T, TOPO >::offsetBy(), posCalculator_, and recHits_.

Referenced by mainSearch().

 {
   
   DetId thisDet;
   std::vector< std::pair<DetId, float> > dets;
   dets.clear();
   EcalRecHitCollection::const_iterator hit;
   
   for (int dx = -2; dx < 3; ++dx)
     {
       for (int dy = -2; dy < 3; ++ dy)
     {
       thisDet = navigator.offsetBy(dx, dy);
       navigator.home();
       
       if (thisDet != DetId(0))
         {
           hit = recHits_->find(thisDet);
           if (hit != recHits_->end()) 
         {
           dets.push_back( std::pair<DetId, float>(thisDet, 1.) ); // by default hit energy fraction is set to 1
         }
         }
     }
     }
   
   // compute coefficient to turn E into Et 
   Point pos = posCalculator_.Calculate_Location(dets, hits, geometry);
   return  1/cosh(pos.eta());
 
 }

double HybridClusterAlgo::et25	(	EcalBarrelNavigatorHT &	navigator,
		const EcalRecHitCollection *	hits,
		const CaloSubdetectorGeometry *	geometry
	)

private

Definition at line 626 of file HybridClusterAlgo.cc.

References PositionCalc::Calculate_Location(), edm::SortedCollection< T, SORT >::end(), CastorDataFrameFilter_impl::energySum(), edm::SortedCollection< T, SORT >::find(), CaloNavigator< T, TOPO >::home(), CaloNavigator< T, TOPO >::offsetBy(), posCalculator_, and recHits_.

Referenced by mainSearch().

 {
 
     DetId thisDet;
     std::vector< std::pair<DetId, float> > dets;
     dets.clear();
     EcalRecHitCollection::const_iterator hit;
     double energySum = 0.0;
 
     for (int dx = -2; dx < 3; ++dx)
     {
         for (int dy = -2; dy < 3; ++ dy)
         {
             //std::cout << "dx, dy " << dx << ", " << dy << std::endl;
             thisDet = navigator.offsetBy(dx, dy);
             navigator.home();
 
             if (thisDet != DetId(0))
             {
                 hit = recHits_->find(thisDet);
                 if (hit != recHits_->end()) 
                 {
                     dets.push_back( std::pair<DetId, float>(thisDet, 1.) ); // by default hit energy fraction is set to 1
                     energySum += hit->energy();
                 }
             }
         }
     }
 
     // convert it to ET
     //std::cout << "dets.size(), energySum: " << dets.size() << ", " << energySum << std::endl;
     Point pos = posCalculator_.Calculate_Location(dets, hits, geometry);
     double et = energySum/cosh(pos.eta());
     return et;
 
 }

void HybridClusterAlgo::mainSearch	(	const EcalRecHitCollection *	hits,
		const CaloSubdetectorGeometry *	geometry
	)

Definition at line 181 of file HybridClusterAlgo.cc.

References BremRecoveryPhiRoadAlgo::barrelPhiRoad(), PositionCalc::Calculate_Location(), clustered_, reco::CaloID::DET_ECAL_BARREL, dynamicEThres_, dynamicPhiRoad_, e2Et(), Eseed, et25(), eThres_, eThresA_, eThresB_, CaloNavigator< T, TOPO >::home(), reco::CaloCluster::hybrid, i, j, relval_steps::k, LogTrace, makeDomino(), HLT_25ns14e33_v1_cff::navigator, CaloNavigator< T, TOPO >::north(), DetId::null(), phiRoadAlgo_, phiSteps_, posCalculator_, HLT_25ns14e33_v3_cff::recHits, seedClus_, seeds, CaloNavigator< T, TOPO >::south(), mathSSE::sqrt(), groupFilesInBlocks::temp, topo_, useddetids, UseEtForXi, and Xi.

Referenced by makeClusters().

 {
   LogTrace("EcalClusters") << "HybridClusterAlgo Algorithm - looking for clusters" ;
   LogTrace("EcalClusters") << "Found the following clusters:" ;
     
     // Loop over seeds:
     std::vector<EcalRecHit>::iterator it;
     int clustercounter=0;
 
     for (it = seeds.begin(); it != seeds.end(); it++){
         std::vector <reco::BasicCluster> thisseedClusters;
         DetId itID = it->id();
 
         // make sure the current seed has not been used/will not be used in the future:
         std::set<DetId>::iterator seed_in_rechits_it = useddetids.find(itID);
 
         if (seed_in_rechits_it != useddetids.end()) continue;
         //If this seed is already used, then don't use it again.
 
         // output some info on the hit:
         LogTrace("EcalClusters") << "*****************************************************" << "Seed of energy E = " << it->energy() << " @ " << EBDetId(itID) << "*****************************************************" ;
         
 
         //Make a navigator, and set it to the seed cell.
         EcalBarrelNavigatorHT navigator(itID, topo_);
 
         //Now use the navigator to start building dominoes.
 
         //Walking positive in phi:
         std::vector <double> dominoEnergyPhiPlus;  //Here I will store the results of the domino sums
         std::vector <std::vector <EcalRecHit> > dominoCellsPhiPlus; //These are the actual EcalRecHit for dominos.
 
         //Walking negative in phi
         std::vector <double> dominoEnergyPhiMinus;  //Here I will store the results of the domino sums
         std::vector <std::vector <EcalRecHit> > dominoCellsPhiMinus; //These are the actual EcalRecHit for dominos.
 
         //The two sets together.
         std::vector <double> dominoEnergy;  //Here I will store the results of the domino sums
         std::vector <std::vector <EcalRecHit> > dominoCells; //These are the actual EcalRecHit for dominos.
 
         //First, the domino about the seed:
         std::vector <EcalRecHit> initialdomino;
         double e_init = makeDomino(navigator, initialdomino);
         if (e_init < Eseed) continue;
 
     LogTrace("EcalClusters") << "Made initial domino" ;
     
         //
         // compute the phi road length
         double phiSteps;
         double et5x5 = et25(navigator, hits, geometry);
         if (dynamicPhiRoad_ && e_init > 0)
         {
             phiSteps = phiRoadAlgo_->barrelPhiRoad(et5x5);
             navigator.home();
         } else phiSteps = phiSteps_;
 
         //Positive phi steps.
         for (int i=0;i<phiSteps;++i){
             //remember, this always increments the current position of the navigator.
             DetId centerD = navigator.north();
             if (centerD.null())
                 continue;
                         LogTrace("EcalClusters") << "Step ++" << i << " @ " << EBDetId(centerD) ;
             EcalBarrelNavigatorHT dominoNav(centerD, topo_);
 
             //Go get the new domino.
             std::vector <EcalRecHit> dcells;
             double etemp = makeDomino(dominoNav, dcells);
 
             //save this information
             dominoEnergyPhiPlus.push_back(etemp);
             dominoCellsPhiPlus.push_back(dcells);
         }
         LogTrace("EcalClusters") << "Got positive dominos" ;
         //return to initial position
         navigator.home();
 
         //Negative phi steps.
         for (int i=0;i<phiSteps;++i){
             //remember, this always decrements the current position of the navigator.
             DetId centerD = navigator.south();
             if (centerD.null())
                 continue;
 
             LogTrace("EcalClusters") << "Step --" << i << " @ " << EBDetId(centerD) ;
             EcalBarrelNavigatorHT dominoNav(centerD, topo_);
 
             //Go get the new domino.
             std::vector <EcalRecHit> dcells;
             double etemp = makeDomino(dominoNav, dcells);
 
             //save this information
             dominoEnergyPhiMinus.push_back(etemp);
             dominoCellsPhiMinus.push_back(dcells);
         }
         
         LogTrace("EcalClusters") << "Got negative dominos: " ;
     
         //Assemble this information:
         for (int i=int(dominoEnergyPhiMinus.size())-1;i >= 0;--i){
             dominoEnergy.push_back(dominoEnergyPhiMinus[i]);
             dominoCells.push_back(dominoCellsPhiMinus[i]);
         }
         dominoEnergy.push_back(e_init);
         dominoCells.push_back(initialdomino);
         for (int i=0;i<int(dominoEnergyPhiPlus.size());++i){
             dominoEnergy.push_back(dominoEnergyPhiPlus[i]);
             dominoCells.push_back(dominoCellsPhiPlus[i]);
         }
 
         //Ok, now I have all I need in order to go ahead and make clusters.
         LogTrace("EcalClusters") << "Dumping domino energies: " ;
         
 
         //for (int i=0;i<int(dominoEnergy.size());++i){
         //       LogTrace("EcalClusters") << "Domino: " << i << " E: " << dominoEnergy[i] ;
         //      }
         //Identify the peaks in this set of dominos:
         //Peak==a domino whose energy is greater than the two adjacent dominos.
         //thus a peak in the local sense.
         double etToE(1.);
         if(!UseEtForXi){
           etToE = 1./e2Et(navigator, hits, geometry);
         }
         std::vector <int> PeakIndex;
         for (int i=1;i<int(dominoEnergy.size())-1;++i){
             if (dominoEnergy[i] > dominoEnergy[i-1]
                     && dominoEnergy[i] >= dominoEnergy[i+1]
                 && dominoEnergy[i] > sqrt( Eseed*Eseed + Xi*Xi*et5x5*et5x5*etToE*etToE) )
               { // Eseed increases ~proportiaonally to et5x5
                 PeakIndex.push_back(i);
               }
         }
 
         LogTrace("EcalClusters") << "Found: " << PeakIndex.size() << " peaks." ;
     
         //Order these peaks by energy:
         for (int i=0;i<int(PeakIndex.size());++i){
             for (int j=0;j<int(PeakIndex.size())-1;++j){
                 if (dominoEnergy[PeakIndex[j]] < dominoEnergy[PeakIndex[j+1]]){
                     int ihold = PeakIndex[j+1];
                     PeakIndex[j+1] = PeakIndex[j];
                     PeakIndex[j] = ihold;
                 }
             }
         }
 
         std::vector<int> OwnerShip;
         std::vector<double> LumpEnergy;
         for (int i=0;i<int(dominoEnergy.size());++i) OwnerShip.push_back(-1);
 
         //Loop over peaks.  
         double eThres = eThres_;
         double e5x5 = 0.0;
         for (int i = 0; i < int(PeakIndex.size()); ++i)
         {
 
             int idxPeak = PeakIndex[i];
             OwnerShip[idxPeak] = i;
             double lump = dominoEnergy[idxPeak];
 
             // compute eThres for this peak
             // if set to dynamic (otherwise uncanged from
             // fixed setting
             if (dynamicEThres_) {
 
                 //std::cout << "i : " << i << " idxPeak " << idxPeak << std::endl;
                 //std::cout << "    the dominoEnergy.size() = " << dominoEnergy.size() << std::endl;
                 // compute e5x5 for this seed crystal
                 //std::cout << "idxPeak, phiSteps " << idxPeak << ", " << phiSteps << std::endl;
                 e5x5 = lump;
                 //std::cout << "lump " << e5x5 << std::endl;
                 if ((idxPeak + 1) < (int)dominoEnergy.size()) e5x5 += dominoEnergy[idxPeak + 1];
                 //std::cout << "+1 " << e5x5 << std::endl;
                 if ((idxPeak + 2) < (int)dominoEnergy.size()) e5x5 += dominoEnergy[idxPeak + 2];
                 //std::cout << "+2 " << e5x5 << std::endl;
                 if ((idxPeak - 1) > 0) e5x5 += dominoEnergy[idxPeak - 1];
                 //std::cout << "-1 " << e5x5 << std::endl;
                 if ((idxPeak - 2) > 0) e5x5 += dominoEnergy[idxPeak - 2];
                 //std::cout << "-2 " << e5x5 << std::endl;
                 // compute eThres
                 eThres = (eThresA_ * e5x5) + eThresB_;   
                 //std::cout << eThres << std::endl;
                 //std::cout << std::endl;
             }
 
             //Loop over adjacent dominos at higher phi
             for (int j=idxPeak+1;j<int(dominoEnergy.size());++j){
                 if (OwnerShip[j]==-1 && 
                         dominoEnergy[j] > eThres
                         && dominoEnergy[j] <= dominoEnergy[j-1]){
                     OwnerShip[j]= i;
                     lump+=dominoEnergy[j];
                 }
                 else{
                     break;
                 }
             }
             //loop over adjacent dominos at lower phi.  Sum up energy of lumps.
             for (int j=idxPeak-1;j>=0;--j){
                 if (OwnerShip[j]==-1 && 
                         dominoEnergy[j] > eThres
                         && dominoEnergy[j] <= dominoEnergy[j+1]){
                     OwnerShip[j]= i;
                     lump+=dominoEnergy[j];
                 }
                 else{
                     break;
                 }
             }
             LumpEnergy.push_back(lump);
         }
 
         //Make the basic clusters:
         for (int i=0;i<int(PeakIndex.size());++i){
             bool HasSeedCrystal = false;
             //One cluster for each peak.
             std::vector<EcalRecHit> recHits;
             std::vector< std::pair<DetId, float> > dets;
             int nhits=0;
             for (int j=0;j<int(dominoEnergy.size());++j){   
                 if (OwnerShip[j] == i){
                     std::vector <EcalRecHit> temp = dominoCells[j];
                     for (int k=0;k<int(temp.size());++k){
                         dets.push_back( std::pair<DetId, float>(temp[k].id(), 1.) ); // by default energy fractions are 1
                         if (temp[k].id()==itID)
                             HasSeedCrystal = true;
                         recHits.push_back(temp[k]);
                         nhits++;
                     }
                 }  
             }
                         LogTrace("EcalClusters") << "Adding a cluster with: " << nhits;
             LogTrace("EcalClusters") << "total E: " << LumpEnergy[i];
             LogTrace("EcalClusters") << "total dets: " << dets.size() ;
     
             //Get Calorimeter position
             Point pos = posCalculator_.Calculate_Location(dets,hits,geometry);
 
             //double totChi2=0;
             //double totE=0;
             std::vector<std::pair<DetId, float> > usedHits;
             for (int blarg=0;blarg<int(recHits.size());++blarg){
                 //totChi2 +=0;
                 //totE+=recHits[blarg].energy();
                 usedHits.push_back(std::make_pair<DetId, float>(recHits[blarg].id(), 1.0));
                 useddetids.insert(recHits[blarg].id());
             }
 
             //if (totE>0)
             //totChi2/=totE;
 
             if (HasSeedCrystal) {
                 // note that this "basiccluster" has the seed crystal of the hyrbid, so record it
                 seedClus_.push_back(reco::BasicCluster(LumpEnergy[i], pos, 
                     reco::CaloID(reco::CaloID::DET_ECAL_BARREL), usedHits, 
                     reco::CaloCluster::hybrid, itID));
                 // and also add to the vector of clusters that will be used in constructing
                 // the supercluster
                                 thisseedClusters.push_back(reco::BasicCluster(LumpEnergy[i], pos,                                                               reco::CaloID(reco::CaloID::DET_ECAL_BARREL), 
                                      usedHits, reco::CaloCluster::hybrid, itID));
             }
             else {
                 // note that if this "basiccluster" is not the one that seeded the hybrid, 
                 // the seed crystal is unset in this entry in the vector of clusters that will
                 // be used in constructing the super cluster
                 thisseedClusters.push_back(reco::BasicCluster(LumpEnergy[i], pos,                                                   reco::CaloID(reco::CaloID::DET_ECAL_BARREL), 
                      usedHits, reco::CaloCluster::hybrid));
             }
         }
 
 
         // Make association so that superclusters can be made later.
         // but only if some BasicClusters have been found...
         if (thisseedClusters.size() > 0) 
         {
             clustered_.insert(std::make_pair(clustercounter, thisseedClusters));
             clustercounter++;
         }
     }//Seed loop
 
 }// end of mainSearch

void HybridClusterAlgo::makeClusters	(	const EcalRecHitCollection *	recColl,
		const CaloSubdetectorGeometry *	geometry,
		reco::BasicClusterCollection &	basicClusters,
		const EcalSeverityLevelAlgo *	sevLv,
		bool	regional = `false`,
		const std::vector< EcalEtaPhiRegion > &	regions = `std::vector<EcalEtaPhiRegion>()`
	)

Definition at line 57 of file HybridClusterAlgo.cc.

References edm::SortedCollection< T, SORT >::begin(), clustered_, eb_st, edm::SortedCollection< T, SORT >::end(), ET, excludedCrys_, excludeFromCluster_, spr::find(), CaloCellGeometry::getPosition(), j, EcalRecHit::kGood, LogTrace, mainSearch(), position, recHits_, HLT_25ns14e33_v3_cff::region, seedClus_, seeds, EcalSeverityLevelAlgo::severityLevel(), funct::sin(), python.multivaluedict::sort(), PV3DBase< T, PVType, FrameType >::theta(), useddetids, v_chstatus_, and v_severitylevel_.

Referenced by EgammaHLTHybridClusterProducer::produce(), and HybridClusterProducer::produce().

 {
   //clear vector of seeds
   seeds.clear();
   //clear flagged crystals
   excludedCrys_.clear();
   //clear map of supercluster/basiccluster association
   clustered_.clear();
   //clear set of used detids
   useddetids.clear();
   //clear vector of seed clusters
   seedClus_.clear();
   //Pass in a pointer to the collection.
   recHits_ = recColl;
   
   LogTrace("EcalClusters") << "Cleared vectors, starting clusterization..." ;
   LogTrace("EcalClusters") << " Purple monkey aardvark." ;
   //
   
   int nregions=0;
   if(regional) nregions=regions.size();
   
   if(!regional || nregions) {
     
     EcalRecHitCollection::const_iterator it;
     
     for (it = recHits_->begin(); it != recHits_->end(); it++){
       
       //Make the vector of seeds that we're going to use.
       //One of the few places position is used, needed for ET calculation.    
       const CaloCellGeometry *this_cell = (*geometry).getGeometry(it->id());
       GlobalPoint position = this_cell->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) {
 
     //Must pass seed threshold
     float ET = it->energy() * sin(position.theta());
 
     LogTrace("EcalClusters") << "Seed crystal: " ;
 
     if (ET > eb_st) {
       
 
       // avoid seeding for anomalous channels   
       if(! it->checkFlag(EcalRecHit::kGood)){ // if rechit is good, no need for further checks
         if (it->checkFlags( v_chstatus_ )) {        
           if (excludeFromCluster_) excludedCrys_.insert(it->id());
           continue; // the recHit has to be excluded from seeding       
         }       
       }
 
 
       int severityFlag =  sevLv->severityLevel( it->id(), *recHits_);
       std::vector<int>::const_iterator sit = std::find( v_severitylevel_.begin(), v_severitylevel_.end(), severityFlag);
       LogTrace("EcalClusters") << "found flag: " << severityFlag ;
        
       
       if (sit!=v_severitylevel_.end()){ 
         if (excludeFromCluster_)
           excludedCrys_.insert(it->id());
         continue;
       }
       seeds.push_back(*it);
      
       LogTrace("EcalClusters") << "Seed ET: " << ET ;
       LogTrace("EcalClsuters") << "Seed E: " << it->energy() ;
     }
       }
     }
     
   }
   
   
   //Yay sorting.
   LogTrace("EcalClusters") << "Built vector of seeds, about to sort them..." ;
   
   //Needs three argument sort with seed comparison operator
   sort(seeds.begin(), seeds.end(), less_mag());
   
   LogTrace("EcalClusters") << "done" ;
   
   //Now to do the work.
   LogTrace("EcalClusters") << "About to call mainSearch...";
   mainSearch(recColl,geometry);
   LogTrace("EcalClusters") << "done" ;
   
   //Hand the basicclusters back to the producer.  It has to 
   //put them in the event.  Then we can make superclusters.
   std::map<int, reco::BasicClusterCollection>::iterator bic; 
   for (bic= clustered_.begin();bic!=clustered_.end();bic++){
     reco::BasicClusterCollection bl = bic->second;
     for (int j=0;j<int(bl.size());++j){
       basicClusters.push_back(bl[j]);
     }
   }
   
   //Yay more sorting.
   sort(basicClusters.rbegin(), basicClusters.rend(), ClusterEtLess() );
   //Done!
   LogTrace("EcalClusters") << "returning to producer. ";
   
 }

double HybridClusterAlgo::makeDomino	(	EcalBarrelNavigatorHT &	navigator,
		std::vector< EcalRecHit > &	cells
	)

Definition at line 542 of file HybridClusterAlgo.cc.

References CaloNavigator< T, TOPO >::east(), edm::SortedCollection< T, SORT >::end(), EcalRecHit::energy(), Ewing, excludedCrys_, edm::SortedCollection< T, SORT >::find(), CaloNavigator< T, TOPO >::home(), CaloNavigator< T, TOPO >::pos(), recHits_, useddetids, and CaloNavigator< T, TOPO >::west().

Referenced by mainSearch().

 {
     //At the beginning of this function, the navigator starts at the middle of the domino,
     //and that's where EcalBarrelNavigator::home() should send it.
     //Walk one crystal in eta to either side of the initial point.  Sum the three cell energies.
     //If the resultant energy is larger than Ewing, then go an additional cell to either side.
     //Returns:  Total domino energy.  Also, stores the cells used to create domino in the vector.
     cells.clear();
     double Etot = 0;
 
     //Ready?  Get the starting cell.
     DetId center = navigator.pos();
     EcalRecHitCollection::const_iterator center_it = recHits_->find(center);
 
     if (center_it!=recHits_->end()){
         EcalRecHit SeedHit = *center_it;
         if (useddetids.find(center) == useddetids.end() && excludedCrys_.find(center)==excludedCrys_.end()){ 
             Etot += SeedHit.energy();
             cells.push_back(SeedHit);
         }
     }
     //One step upwards in Ieta:
     DetId ieta1 = navigator.west();
     EcalRecHitCollection::const_iterator eta1_it = recHits_->find(ieta1);
     if (eta1_it !=recHits_->end()){
         EcalRecHit UpEta = *eta1_it;
         if (useddetids.find(ieta1) == useddetids.end() && excludedCrys_.find(ieta1)==excludedCrys_.end()){
             Etot+=UpEta.energy();
             cells.push_back(UpEta);
         }
     }
 
     //Go back to the middle.
     navigator.home();
 
     //One step downwards in Ieta:
     DetId ieta2 = navigator.east();
     EcalRecHitCollection::const_iterator eta2_it = recHits_->find(ieta2);
     if (eta2_it !=recHits_->end()){
         EcalRecHit DownEta = *eta2_it;
         if (useddetids.find(ieta2)==useddetids.end() && excludedCrys_.find(ieta2)==excludedCrys_.end()){
             Etot+=DownEta.energy();
             cells.push_back(DownEta);
         }
     }
 
     //Now check the energy.  If smaller than Ewing, then we're done.  If greater than Ewing, we have to
     //add two additional cells, the 'wings'
     if (Etot < Ewing) {
         navigator.home(); //Needed even here!!
         return Etot;  //Done!  Not adding 'wings'.
     }
 
     //Add the extra 'wing' cells.  Remember, we haven't sent the navigator home,
     //we're still on the DownEta cell.
     if (ieta2 != DetId(0)){
         DetId ieta3 = navigator.east(); //Take another step downward.
         EcalRecHitCollection::const_iterator eta3_it = recHits_->find(ieta3);
         if (eta3_it != recHits_->end()){
             EcalRecHit DownEta2 = *eta3_it;
             if (useddetids.find(ieta3)==useddetids.end() && excludedCrys_.find(ieta3)==excludedCrys_.end()){
                 Etot+=DownEta2.energy();
                 cells.push_back(DownEta2);
             }
         }
     }
     //Now send the navigator home.
     navigator.home();
     navigator.west(); //Now you're on eta1_it
     if (ieta1 != DetId(0)){
         DetId ieta4 = navigator.west(); //Take another step upward.
         EcalRecHitCollection::const_iterator eta4_it = recHits_->find(ieta4);
         if (eta4_it != recHits_->end()){
             EcalRecHit UpEta2 = *eta4_it;
             if (useddetids.find(ieta4) == useddetids.end() && excludedCrys_.find(ieta4)==excludedCrys_.end()){
                 Etot+=UpEta2.energy();
                 cells.push_back(UpEta2);
             }
         }
     }
     navigator.home();
     return Etot;
 }

reco::SuperClusterCollection HybridClusterAlgo::makeSuperClusters ( const reco::CaloClusterPtrVector & clustersCollection )

Definition at line 465 of file HybridClusterAlgo.cc.

References clustered_, reco::CaloCluster::energy(), i, j, LogTrace, edm::PtrVector< T >::push_back(), fileCollector::seed, seedClus_, edm::PtrVectorBase::size(), and python.multivaluedict::sort().

Referenced by EgammaHLTHybridClusterProducer::produce(), and HybridClusterProducer::produce().

 {
     //Here's what we'll return.
     reco::SuperClusterCollection SCcoll;
 
     //Here's our map iterator that gives us the appropriate association.
     std::map<int, reco::BasicClusterCollection>::iterator mapit;
     for (mapit = clustered_.begin();mapit!=clustered_.end();mapit++){
 
         reco::CaloClusterPtrVector thissc;
         reco::CaloClusterPtr seed;//This is not really a seed, but I need to tell SuperCluster something.
         //So I choose the highest energy basiccluster in the SuperCluster.
 
         std::vector <reco::BasicCluster> thiscoll = mapit->second; //This is the set of BasicClusters in this
         //SuperCluster
 
         double ClusterE = 0; //Sum of cluster energies for supercluster.
         //Holders for position of this supercluster.
         double posX=0;
         double posY=0;
         double posZ=0;
 
         //Loop over this set of basic clusters, find their references, and add them to the
         //supercluster.  This could be somehow more efficient.
 
         for (int i=0;i<int(thiscoll.size());++i){
             reco::BasicCluster thisclus = thiscoll[i]; //The Cluster in question.
             for (int j=0;j<int(clustersCollection.size());++j){
                 //Find the appropriate cluster from the list of references
                 reco::BasicCluster cluster_p = *clustersCollection[j];
                 if (thisclus== cluster_p){ //Comparison based on energy right now.
                     thissc.push_back(clustersCollection[j]);
                     bool isSeed = false;
                     for (int qu=0;qu<int(seedClus_.size());++qu){
                         if (cluster_p == seedClus_[qu])
                             isSeed = true;
                     }
                     if (isSeed) seed = clustersCollection[j];
 
                     ClusterE += cluster_p.energy();
                     posX += cluster_p.energy() * cluster_p.position().X();
                     posY += cluster_p.energy() * cluster_p.position().Y();
                     posZ += cluster_p.energy() * cluster_p.position().Z();
 
                 }
             }//End loop over finding references.
         }//End loop over clusters.
 
         posX /= ClusterE;
         posY /= ClusterE;
         posZ /= ClusterE;
 
         /* //This part is moved to EgammaSCEnergyCorrectionAlgo
            double preshowerE = 0.;
            double phiWidth = 0.;
            double etaWidth = 0.;
         //Calculate phiWidth & etaWidth for SuperClusters
         reco::SuperCluster suCl(ClusterE, math::XYZPoint(posX, posY, posZ), seed, thissc, preshowerE, phiWidth, etaWidth);
         SCShape_->Calculate_Covariances(suCl);
         phiWidth = SCShape_->phiWidth();
         etaWidth = SCShape_->etaWidth();
         //Assign phiWidth & etaWidth to SuperCluster as data members
         suCl.setPhiWidth(phiWidth);
         suCl.setEtaWidth(etaWidth);
          */
 
         reco::SuperCluster suCl(ClusterE, math::XYZPoint(posX, posY, posZ), seed, thissc);
         SCcoll.push_back(suCl);
 
                 LogTrace("EcalClusters") << "Super cluster sum: " << ClusterE ;
                 LogTrace("EcalClusters") << "Made supercluster with energy E: " << suCl.energy() ;
         
     }//end loop over map
     sort(SCcoll.rbegin(), SCcoll.rend(), ClusterEtLess());
     return SCcoll;
 }

void HybridClusterAlgo::setDynamicPhiRoad ( const edm::ParameterSet & bremRecoveryPset )

inline

Definition at line 154 of file HybridClusterAlgo.h.

References dynamicPhiRoad_, and phiRoadAlgo_.

   {
      dynamicPhiRoad_ = true;
      phiRoadAlgo_ = new BremRecoveryPhiRoadAlgo(bremRecoveryPset);
   }

Member Data Documentation

std::map<int, std::vector<reco::BasicCluster> > HybridClusterAlgo::clustered_

private

Definition at line 104 of file HybridClusterAlgo.h.

Referenced by mainSearch(), makeClusters(), and makeSuperClusters().

bool HybridClusterAlgo::dynamicEThres_

private

Definition at line 80 of file HybridClusterAlgo.h.

Referenced by HybridClusterAlgo(), and mainSearch().

bool HybridClusterAlgo::dynamicPhiRoad_

private

Definition at line 77 of file HybridClusterAlgo.h.

Referenced by HybridClusterAlgo(), mainSearch(), setDynamicPhiRoad(), and ~HybridClusterAlgo().

double HybridClusterAlgo::eb_st

private

Definition at line 40 of file HybridClusterAlgo.h.

Referenced by makeClusters().

double HybridClusterAlgo::Eseed

private

Definition at line 65 of file HybridClusterAlgo.h.

Referenced by HybridClusterAlgo(), and mainSearch().

double HybridClusterAlgo::eThres_

private

Definition at line 60 of file HybridClusterAlgo.h.

Referenced by mainSearch().

double HybridClusterAlgo::eThresA_

private

Definition at line 61 of file HybridClusterAlgo.h.

Referenced by HybridClusterAlgo(), and mainSearch().

double HybridClusterAlgo::eThresB_

private

Definition at line 62 of file HybridClusterAlgo.h.

Referenced by HybridClusterAlgo(), and mainSearch().

double HybridClusterAlgo::Ewing

private

Definition at line 74 of file HybridClusterAlgo.h.

Referenced by makeDomino().

std::set<DetId> HybridClusterAlgo::excludedCrys_

private

Definition at line 118 of file HybridClusterAlgo.h.

Referenced by makeClusters(), and makeDomino().

bool HybridClusterAlgo::excludeFromCluster_

private

Definition at line 117 of file HybridClusterAlgo.h.

Referenced by makeClusters().

BremRecoveryPhiRoadAlgo* HybridClusterAlgo::phiRoadAlgo_

private

Definition at line 57 of file HybridClusterAlgo.h.

Referenced by mainSearch(), setDynamicPhiRoad(), and ~HybridClusterAlgo().

int HybridClusterAlgo::phiSteps_

private

Definition at line 45 of file HybridClusterAlgo.h.

Referenced by mainSearch().

PositionCalc HybridClusterAlgo::posCalculator_

private

Definition at line 107 of file HybridClusterAlgo.h.

Referenced by e2Et(), et25(), HybridClusterAlgo(), and mainSearch().

const EcalRecHitCollection* HybridClusterAlgo::recHits_

private

Definition at line 87 of file HybridClusterAlgo.h.

Referenced by e2Et(), et25(), makeClusters(), and makeDomino().

std::vector<reco::BasicCluster> HybridClusterAlgo::seedClus_

private

Definition at line 101 of file HybridClusterAlgo.h.

Referenced by mainSearch(), makeClusters(), and makeSuperClusters().

std::vector<EcalRecHit> HybridClusterAlgo::seeds

private

Definition at line 98 of file HybridClusterAlgo.h.

Referenced by mainSearch(), and makeClusters().

float HybridClusterAlgo::severityRecHitThreshold_

private

Definition at line 114 of file HybridClusterAlgo.h.

float HybridClusterAlgo::severitySpikeThreshold_

private

Definition at line 115 of file HybridClusterAlgo.h.

EcalBarrelHardcodedTopology* HybridClusterAlgo::topo_

private

Definition at line 90 of file HybridClusterAlgo.h.

Referenced by HybridClusterAlgo(), mainSearch(), and ~HybridClusterAlgo().

std::set<DetId> HybridClusterAlgo::useddetids

private

Definition at line 95 of file HybridClusterAlgo.h.

Referenced by mainSearch(), makeClusters(), and makeDomino().

bool HybridClusterAlgo::UseEtForXi

private

Definition at line 71 of file HybridClusterAlgo.h.

Referenced by HybridClusterAlgo(), and mainSearch().

std::vector<int> HybridClusterAlgo::v_chstatus_

private

Definition at line 110 of file HybridClusterAlgo.h.

Referenced by HybridClusterAlgo(), and makeClusters().

std::vector<int> HybridClusterAlgo::v_severitylevel_

private

Definition at line 113 of file HybridClusterAlgo.h.

Referenced by makeClusters().

double HybridClusterAlgo::Xi

private

Definition at line 68 of file HybridClusterAlgo.h.

Referenced by HybridClusterAlgo(), and mainSearch().

Public Member Functions

Private Types

Private Member Functions

Private Attributes

Detailed Description

Member Typedef Documentation

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation