CMSSW Reference Manual

typedef std::set<DetId> PreshowerClusterAlgo::HitsID

Definition at line 34 of file PreshowerClusterAlgo.h.

typedef math::XYZPoint PreshowerClusterAlgo::Point

Definition at line 31 of file PreshowerClusterAlgo.h.

typedef std::map<DetId, EcalRecHit> PreshowerClusterAlgo::RecHitsMap

Definition at line 33 of file PreshowerClusterAlgo.h.

Enumerator:

pDEBUG
pINFO
pERROR

Definition at line 29 of file PreshowerClusterAlgo.h.

00029 { pDEBUG = 0, pINFO = 1, pERROR = 2 };

PreshowerClusterAlgo::PreshowerClusterAlgo ( ) [inline]

Definition at line 36 of file PreshowerClusterAlgo.h.

00036                           : 
00037    preshStripEnergyCut_(0.), preshClusterEnergyCut_(0.), preshSeededNstr_(15), debugLevel_(pINFO)
00038    {}

PreshowerClusterAlgo::PreshowerClusterAlgo	(	double	stripEnergyCut,
		double	clusterEnergyCut,
		int	nStripCut,
		DebugLevel	debugLevel = `pINFO`
	)			`[inline]`

Definition at line 40 of file PreshowerClusterAlgo.h.

00040                                                                                                                       :
00041    preshStripEnergyCut_(stripEnergyCut), preshClusterEnergyCut_(clusterEnergyCut), preshSeededNstr_(nStripCut), debugLevel_(debugLevel)
00042    {}

PreshowerClusterAlgo::~PreshowerClusterAlgo ( ) [inline]

Definition at line 44 of file PreshowerClusterAlgo.h.

00044 {};

void PreshowerClusterAlgo::findRoad	(	ESDetId	strip,
		EcalPreshowerNavigator	theESNav,
		int	plane
	)

Definition at line 267 of file PreshowerClusterAlgo.cc.

References GenMuonPlsPt100GeV_cfg::cout, debugLevel_, CaloNavigator< T >::east(), lat::endl(), CaloNavigator< T >::home(), CaloNavigator< T >::north(), pDEBUG, pINFO, preshSeededNstr_, road_2d, CaloNavigator< T >::setHome(), CaloNavigator< T >::south(), and CaloNavigator< T >::west().

Referenced by makeOneCluster().

00267                                                                                              {
00268   
00269   if ( strip == ESDetId(0) ) return;
00270 
00271    ESDetId next;
00272    theESNav.setHome(strip);
00273 // First, add a central strip to the road 
00274    road_2d.push_back(strip);   
00275    
00276    if ( debugLevel_ <= pINFO ) std::cout << "findRoad starts from strip " << strip << std::endl;  
00277    if (plane == 1) {
00278      // east road
00279      int n_east= 0;
00280      if ( debugLevel_ == pDEBUG ) std::cout << " Go to the East " <<  std::endl;   
00281      while ( ((next=theESNav.east()) != ESDetId(0) && next != strip) ) {
00282         if ( debugLevel_ == pDEBUG ) std::cout << "East: " << n_east << " current strip is " << next << std::endl;  
00283         road_2d.push_back(next);   
00284         ++n_east;  
00285         if (n_east == preshSeededNstr_) break; 
00286      }
00287      // west road
00288      int n_west= 0;
00289      if ( debugLevel_ == pDEBUG ) std::cout << " Go to the West " <<  std::endl;
00290      theESNav.home();
00291      while ( ((next=theESNav.west()) != ESDetId(0) && next != strip )) {
00292         if ( debugLevel_ == pDEBUG ) std::cout << "West: " << n_west << " current strip is " << next << std::endl;  
00293         road_2d.push_back(next);   
00294         ++n_west;  
00295         if (n_west == preshSeededNstr_) break; 
00296      }
00297      if ( debugLevel_ == pDEBUG ) std::cout << "Total number of strips found in the road at 1-st plane is " << n_east+n_west << std::endl;
00298   } 
00299   else if (plane == 2) {
00300     // north road
00301     int n_north= 0;
00302     if ( debugLevel_ == pDEBUG ) std::cout << " Go to the North " <<  std::endl;
00303     while ( ((next=theESNav.north()) != ESDetId(0) && next != strip) ) {       
00304        if ( debugLevel_ == pDEBUG ) std::cout << "North: " << n_north << " current strip is " << next << std::endl; 
00305        road_2d.push_back(next);   
00306        ++n_north;  
00307        if (n_north == preshSeededNstr_) break; 
00308     }
00309     // south road
00310     int n_south= 0;
00311     if ( debugLevel_ == pDEBUG ) std::cout << " Go to the South " <<  std::endl;
00312     theESNav.home();
00313     while ( ((next=theESNav.south()) != ESDetId(0) && next != strip) ) {
00314        if ( debugLevel_ == pDEBUG ) std::cout << "South: " << n_south << " current strip is " << next << std::endl;      
00315        road_2d.push_back(next);   
00316        ++n_south;  
00317        if (n_south == preshSeededNstr_) break; 
00318     }
00319     if ( debugLevel_ == pDEBUG ) std::cout << "Total number of strips found in the road at 2-nd plane is " << n_south+n_north << std::endl;
00320   } 
00321   else {
00322     if ( debugLevel_ == pDEBUG ) std::cout << " Wrong plane number, null cluster will be returned! " << std::endl;    
00323   } // end of if
00324 
00325   theESNav.home();
00326 }

bool PreshowerClusterAlgo::goodStrip ( RecHitsMap::iterator candidate_it )

Definition at line 244 of file PreshowerClusterAlgo.cc.

References GenMuonPlsPt100GeV_cfg::cout, debugLevel_, lat::endl(), pDEBUG, preshStripEnergyCut_, rechits_map, and used_s.

Referenced by makeOneCluster().

00245 {
00246   if ( debugLevel_ == pDEBUG ) {
00247     if ( used_s->find(candidate_it->first) != used_s->end()) 
00248         std::cout << " This strip is in use " << std::endl;    
00249     if (candidate_it == rechits_map->end() )
00250         std::cout << " No such a strip in rechits_map " << std::endl; 
00251     if (candidate_it->second.energy() <= preshStripEnergyCut_)
00252         std::cout << " Strip energy " << candidate_it->second.energy() <<" is below threshold " << std::endl; 
00253   }
00254   // crystal should not be included...
00255   if ( (used_s->find(candidate_it->first) != used_s->end())  ||       // ...if it already belongs to a cluster
00256        (candidate_it == rechits_map->end() )                    ||       // ...if it corresponds to a hit
00257        (candidate_it->second.energy() <= preshStripEnergyCut_ ) )   // ...if it has a negative or zero energy
00258     {
00259       return false;
00260     }
00261 
00262   return true;
00263 }

reco::PreshowerCluster PreshowerClusterAlgo::makeOneCluster	(	ESDetId	strip,
		HitsID *	used_strips,
		RecHitsMap *	rechits_map,
		const CaloSubdetectorGeometry *&	geometry_p,
		CaloSubdetectorTopology *&	topology_p
	)

Definition at line 11 of file PreshowerClusterAlgo.cc.

References edm::SortedCollection< T, SORT >::begin(), GenMuonPlsPt100GeV_cfg::cout, debugLevel_, dummy, CaloNavigator< T >::east(), edm::SortedCollection< T, SORT >::end(), lat::endl(), reco::CaloCluster::energy(), reco::CaloCluster::eta(), findRoad(), CaloSubdetectorGeometry::getGeometry(), CaloCellGeometry::getPosition(), goodStrip(), CaloNavigator< T >::home(), it, reco::PreshowerCluster::nhits(), CaloNavigator< T >::north(), pDEBUG, reco::CaloCluster::phi(), pINFO, ESDetId::plane(), reco::CaloCluster::position(), preshClusterEnergyCut_, edm::SortedCollection< T, SORT >::push_back(), rechits_map, road_2d, CaloNavigator< T >::setHome(), edm::SortedCollection< T, SORT >::size(), CaloNavigator< T >::south(), ESDetId::strip(), used_s, CaloNavigator< T >::west(), PV3DBase< T, PVType, FrameType >::x(), reco::CaloCluster::x(), PV3DBase< T, PVType, FrameType >::y(), reco::CaloCluster::y(), PV3DBase< T, PVType, FrameType >::z(), and reco::CaloCluster::z().

Referenced by PreshowerClusterProducer::produce().

00016 {
00017   road_2d.clear();
00018 
00019   rechits_map = the_rechitsMap_p;
00020 
00021   used_s = used_strips;
00022 
00023   int plane = strip.plane();
00024 
00025   if ( debugLevel_ <= pINFO ) {
00026     std::cout << "Preshower Seeded Algorithm - looking for clusters" << std::endl;
00027     std::cout << "Preshower is intersected at strip " << strip.strip() << ", at plane " << plane << std::endl;
00028   }
00029 
00030   // create null-cluster
00031   std::vector<DetId> dummy;
00032   Point posi(0,0,0);  
00033   if ( debugLevel_ <= pINFO ) std::cout << " Creating a null-cluster" << std::endl;
00034   reco::PreshowerCluster nullcluster=reco::PreshowerCluster(0.,posi,dummy,plane);
00035 
00036   if ( strip == ESDetId(0) ) return nullcluster;   //works in case of no intersected strip found (e.g. in the Barrel)
00037 
00038   // Collection of cluster strips
00039   EcalRecHitCollection clusterRecHits;
00040   // Map of strips for position calculation
00041   RecHitsMap recHits_pos;
00042 
00043   //Make a navigator, and set it to the strip cell.
00044   EcalPreshowerNavigator navigator(strip, topology_p);
00045   navigator.setHome(strip);
00046  //search for neighbours in the central road
00047   findRoad(strip,navigator,plane);
00048   if ( debugLevel_ <= pINFO ) std::cout << "Total number of strips in the central road: " << road_2d.size() << std::endl;
00049 
00050   if ( plane == 1 ) {
00051      ESDetId strip_north = navigator.north();
00052      findRoad(strip_north,navigator,plane);
00053      navigator.home();
00054      ESDetId strip_south = navigator.south();
00055      findRoad(strip_south,navigator,plane);
00056      navigator.home();
00057   }
00058   if ( plane == 2 ) {
00059      ESDetId strip_east = navigator.east();
00060      findRoad(strip_east,navigator,plane);
00061      navigator.home();
00062      ESDetId strip_west = navigator.west();
00063      findRoad(strip_west,navigator,plane);
00064      navigator.home();
00065   }
00066 
00067   if ( debugLevel_ <= pINFO ) std::cout << "Total number of strips in all three roads: " << road_2d.size() << std::endl;
00068 
00069   // Start clustering from strip with max Energy in the road
00070   float E_max = 0.;
00071   bool found = false;
00072   RecHitsMap::iterator max_it;
00073   // Loop over strips:
00074   std::vector<ESDetId>::iterator itID;
00075   for (itID = road_2d.begin(); itID != road_2d.end(); itID++) {
00076     if ( debugLevel_ == pDEBUG ) std::cout << " ID = " << *itID << std::endl;
00077     RecHitsMap::iterator strip_it = rechits_map->find(*itID);   
00078     //if ( strip_it->second.energy() < 0 ) std::cout << "           ##### E = " << strip_it->second.energy() << std::endl;
00079     if(!goodStrip(strip_it)) continue;
00080     if ( debugLevel_ == pDEBUG ) std::cout << " strip is " << ESDetId(strip_it->first) <<"  E = " << strip_it->second.energy() <<std::endl;
00081     float E = strip_it->second.energy();
00082     if ( E > E_max) {
00083        E_max = E;
00084        found = true;
00085        max_it = strip_it;
00086     }
00087   }
00088   // no hottest strip found ==> null cluster will be returned
00089   if ( !found ) return nullcluster;
00090 
00091   // First, save the hottest strip
00092   clusterRecHits.push_back(max_it->second);  
00093   recHits_pos.insert(std::make_pair(max_it->first, max_it->second));
00094   used_s->insert(max_it->first);
00095   if ( debugLevel_ <= pINFO ) {
00096      std::cout << " Central hottest strip " << ESDetId(max_it->first) << " is saved " << std::endl;
00097      std::cout << " with energy E = " <<  E_max << std::endl;    
00098   }
00099 
00100   // Find positions of adjacent strips:
00101   ESDetId next, strip_1, strip_2;
00102   navigator.setHome(max_it->first);
00103   ESDetId startES = max_it->first;
00104 
00105    if (plane == 1) {
00106     // Save two neighbouring strips to the east
00107     int nadjacents_east = 0;
00108     while ( (next=navigator.east()) != ESDetId(0) && next != startES && nadjacents_east < 2 ) {
00109        ++nadjacents_east;
00110        if ( debugLevel_ == pDEBUG ) std::cout << " Adjacent east #" << nadjacents_east <<": "<< next << std::endl;
00111        RecHitsMap::iterator strip_it = rechits_map->find(next);
00112 
00113        if(!goodStrip(strip_it)) continue;
00114        // Save strip for clustering if it exists, not already in use, and satisfies an energy threshold
00115        clusterRecHits.push_back(strip_it->second);       
00116        // save strip for position calculation
00117        if ( nadjacents_east==1 ) strip_1 = next;
00118        used_s->insert(strip_it->first);
00119        if ( debugLevel_ == pDEBUG ) std::cout << " East adjacent strip # " << nadjacents_east << " is saved with energy E = " 
00120                                               << strip_it->second.energy() << std::endl;             
00121     }
00122     // Save two neighbouring strips to the west
00123     navigator.home();
00124     int nadjacents_west = 0;
00125     while ( (next=navigator.west()) != ESDetId(0) && next != startES && nadjacents_west < 2 ) {
00126        ++nadjacents_west;
00127        if ( debugLevel_ == pDEBUG ) std::cout << " Adjacent west #" << nadjacents_west <<": "<< next << std::endl; 
00128        RecHitsMap::iterator strip_it = rechits_map->find(next);
00129        if(!goodStrip(strip_it)) continue;
00130        clusterRecHits.push_back(strip_it->second);
00131        if ( nadjacents_west==1 ) strip_2 = next;
00132        used_s->insert(strip_it->first);
00133        if ( debugLevel_ == pDEBUG ) std::cout << " West adjacent strip # " << nadjacents_west << " is saved with energy E = " 
00134                                              << strip_it->second.energy() << std::endl;           
00135     }
00136   }
00137   else if (plane == 2) {
00138 
00139   // Save two neighbouring strips to the north
00140     int nadjacents_north = 0;
00141     while ( (next=navigator.north()) != ESDetId(0) && next != startES && nadjacents_north < 2 ) {
00142        ++nadjacents_north;
00143        if ( debugLevel_ == pDEBUG ) std::cout << " Adjacent north #" << nadjacents_north <<": "<< next << std::endl;   
00144        RecHitsMap::iterator strip_it = rechits_map->find(next); 
00145        if(!goodStrip(strip_it)) continue;      
00146        clusterRecHits.push_back(strip_it->second);
00147        if ( nadjacents_north==1 ) strip_1 = next;
00148        used_s->insert(strip_it->first);
00149        if ( debugLevel_ == pDEBUG ) std::cout << " North adjacent strip # " << nadjacents_north << " is saved with energy E = " 
00150                                              << strip_it->second.energy() << std::endl;     
00151     }
00152     // Save two neighbouring strips to the south
00153     navigator.home();
00154     int nadjacents_south = 0;
00155     while ( (next=navigator.south()) != ESDetId(0) && next != startES && nadjacents_south < 2 ) {
00156        ++nadjacents_south;
00157        if ( debugLevel_ == pDEBUG ) std::cout << " Adjacent south #" << nadjacents_south <<": "<< next << std::endl;   
00158        RecHitsMap::iterator strip_it = rechits_map->find(next);   
00159        if(!goodStrip(strip_it)) continue;      
00160        clusterRecHits.push_back(strip_it->second);
00161        if ( nadjacents_south==1 ) strip_2 = next;
00162        used_s->insert(strip_it->first);
00163        if ( debugLevel_ == pDEBUG ) std::cout << " South adjacent strip # " << nadjacents_south << " is saved with energy E = " 
00164                                              << strip_it->second.energy() << std::endl;     
00165     }
00166   }
00167   else {
00168     std::cout << " Wrong plane number" << plane <<", null cluster will be returned! " << std::endl;
00169     return nullcluster;
00170   } // end of if
00171   if ( debugLevel_ <=pINFO ) std::cout << " Total size of clusterRecHits is " << clusterRecHits.size() << std::endl;
00172   if ( debugLevel_ <=pINFO ) std::cout << " Two adjacent strips for position calculation are: " 
00173                                       << strip_1 <<" and " << strip_2 << std::endl; 
00174 
00175   // strips for position calculation
00176   RecHitsMap::iterator strip_it1, strip_it2;
00177   if ( strip_1 != ESDetId(0)) {
00178     strip_it1 = rechits_map->find(strip_1);
00179     recHits_pos.insert(std::make_pair(strip_it1->first, strip_it1->second));  
00180   }
00181   if ( strip_2 != ESDetId(0) ) {
00182     strip_it2 = rechits_map->find(strip_2);
00183     recHits_pos.insert(std::make_pair(strip_it2->first, strip_it2->second));  
00184   }
00185 
00186   RecHitsMap::iterator cp;
00187   double energy_pos = 0;
00188   double x_pos = 0;
00189   double y_pos = 0;
00190   double z_pos = 0;
00191   for (cp = recHits_pos.begin(); cp!=recHits_pos.end(); cp++ ) {
00192      double E = cp->second.energy();
00193      energy_pos += E; 
00194      const CaloCellGeometry *thisCell = geometry_p->getGeometry(cp->first);
00195      GlobalPoint position = thisCell->getPosition();
00196      x_pos += E * position.x();
00197      y_pos += E * position.y();
00198      z_pos += E * position.z();     
00199   }
00200   if(energy_pos>0.) {
00201      x_pos /= energy_pos;
00202      y_pos /= energy_pos;
00203      z_pos /= energy_pos;
00204   }
00205   Point pos(x_pos,y_pos,z_pos);
00206   if ( debugLevel_ == pDEBUG ) std::cout << " ES Cluster position = " << "(" << x_pos <<","<< y_pos <<","<< z_pos <<")"<< std::endl;
00207 
00208   EcalRecHitCollection::iterator it;
00209   double Eclust = 0;
00210 
00211   if ( debugLevel_ == pINFO ) std::cout << "The found ES cluster strips: " << std::endl;  
00212   std::vector<DetId> usedHits;
00213   for (it=clusterRecHits.begin(); it != clusterRecHits.end(); it++) {
00214      Eclust += it->energy();
00215      usedHits.push_back(it->id());
00216      if ( debugLevel_ == pINFO ) std::cout << ESDetId(it->id()) <<", E = " << it->energy()<<"; ";
00217   }   
00218   if ( debugLevel_ == pINFO ) std::cout << std::endl;
00219 
00220 
00221   // ES cluster is created from vector clusterRecHits
00222   reco::PreshowerCluster cluster=reco::PreshowerCluster(Eclust,pos,usedHits,plane);
00223 
00224   if ( debugLevel_ <= pINFO ) {
00225      std::cout << " ES Cluster is created with " << std::endl;
00226      std::cout << " energy = " << cluster.energy() << std::endl;
00227      std::cout << " (eta,phi) = " << "("<<cluster.eta()<<", "<<cluster.phi()<<")"<< std::endl;
00228      std::cout << " nhits = " << cluster.nhits() << std::endl;
00229      std::cout << " radius = " << cluster.position().r() << std::endl; 
00230      std::cout << " (x,y,z) = " << "(" << cluster.x() <<", "<< cluster.y() <<", "<< cluster.z()<<")"<< std::endl;     
00231   }
00232  
00233   used_strips = used_s;
00234 
00235   // return the cluster if its energy is greater a threshold
00236   if( cluster.energy() > preshClusterEnergyCut_ ) 
00237      return cluster; 
00238   else
00239      return nullcluster;
00240 
00241 }

int PreshowerClusterAlgo::debugLevel_ [private]

Definition at line 61 of file PreshowerClusterAlgo.h.

Referenced by findRoad(), goodStrip(), and makeOneCluster().

double PreshowerClusterAlgo::preshClusterEnergyCut_ [private]

Definition at line 59 of file PreshowerClusterAlgo.h.

Referenced by makeOneCluster().

int PreshowerClusterAlgo::preshSeededNstr_ [private]

Definition at line 60 of file PreshowerClusterAlgo.h.

Referenced by findRoad().

double PreshowerClusterAlgo::preshStripEnergyCut_ [private]

Definition at line 58 of file PreshowerClusterAlgo.h.

Referenced by goodStrip().

RecHitsMap* PreshowerClusterAlgo::rechits_map [private]

Definition at line 66 of file PreshowerClusterAlgo.h.

Referenced by goodStrip(), and makeOneCluster().

std::vector<ESDetId> PreshowerClusterAlgo::road_2d [private]

Definition at line 63 of file PreshowerClusterAlgo.h.

Referenced by findRoad(), and makeOneCluster().

HitsID* PreshowerClusterAlgo::used_s [private]

Definition at line 69 of file PreshowerClusterAlgo.h.

Referenced by goodStrip(), and makeOneCluster().

PreshowerClusterAlgo Class Reference

Public Types

Public Member Functions

Private Attributes

Detailed Description

Member Typedef Documentation

Member Enumeration Documentation

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation


Public Types
enum	DebugLevel { pDEBUG = 0, pINFO = 1, pERROR = 2 }
typedef std::set< DetId >	HitsID
typedef math::XYZPoint	Point
typedef std::map< DetId, EcalRecHit >	RecHitsMap
Public Member Functions
void	findRoad (ESDetId strip, EcalPreshowerNavigator theESNav, int plane)
bool	goodStrip (RecHitsMap::iterator candidate_it)
reco::PreshowerCluster	makeOneCluster (ESDetId strip, HitsID used_strips, RecHitsMap rechits_map, const CaloSubdetectorGeometry &geometry_p, CaloSubdetectorTopology &topology_p)
	PreshowerClusterAlgo (double stripEnergyCut, double clusterEnergyCut, int nStripCut, DebugLevel debugLevel=pINFO)
	PreshowerClusterAlgo ()
	~PreshowerClusterAlgo ()
Private Attributes
int	debugLevel_
double	preshClusterEnergyCut_
int	preshSeededNstr_
double	preshStripEnergyCut_
RecHitsMap *	rechits_map
std::vector< ESDetId >	road_2d
HitsID *	used_s