|
|
|
#include <PreshowerClusterAlgo.h>
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 |
Definition at line 24 of file PreshowerClusterAlgo.h.
| typedef std::set<DetId> PreshowerClusterAlgo::HitsID |
Definition at line 34 of file PreshowerClusterAlgo.h.
Definition at line 31 of file PreshowerClusterAlgo.h.
| typedef std::map<DetId, EcalRecHit> PreshowerClusterAlgo::RecHitsMap |
Definition at line 33 of file PreshowerClusterAlgo.h.
Definition at line 29 of file PreshowerClusterAlgo.h.
| PreshowerClusterAlgo::PreshowerClusterAlgo | ( | ) | [inline] |
Definition at line 36 of file PreshowerClusterAlgo.h.
: preshStripEnergyCut_(0.), preshClusterEnergyCut_(0.), preshSeededNstr_(15), debugLevel_(pINFO) {}
| PreshowerClusterAlgo::PreshowerClusterAlgo | ( | double | stripEnergyCut, |
| double | clusterEnergyCut, | ||
| int | nStripCut, | ||
| DebugLevel | debugLevel = pINFO |
||
| ) | [inline] |
Definition at line 40 of file PreshowerClusterAlgo.h.
: preshStripEnergyCut_(stripEnergyCut), preshClusterEnergyCut_(clusterEnergyCut), preshSeededNstr_(nStripCut), debugLevel_(debugLevel) {}
| PreshowerClusterAlgo::~PreshowerClusterAlgo | ( | ) | [inline] |
Definition at line 44 of file PreshowerClusterAlgo.h.
{};
| void PreshowerClusterAlgo::findRoad | ( | ESDetId | strip, |
| EcalPreshowerNavigator | theESNav, | ||
| int | plane | ||
| ) |
Definition at line 268 of file PreshowerClusterAlgo.cc.
References gather_cfg::cout, debugLevel_, CaloNavigator< T >::east(), CaloNavigator< T >::home(), CaloNavigator< T >::north(), pDEBUG, pINFO, preshSeededNstr_, road_2d, CaloNavigator< T >::setHome(), CaloNavigator< T >::south(), strip(), and CaloNavigator< T >::west().
Referenced by makeOneCluster().
{
if ( strip == ESDetId(0) ) return;
ESDetId next;
theESNav.setHome(strip);
// First, add a central strip to the road
road_2d.push_back(strip);
if ( debugLevel_ <= pINFO ) std::cout << "findRoad starts from strip " << strip << std::endl;
if (plane == 1) {
// east road
int n_east= 0;
if ( debugLevel_ == pDEBUG ) std::cout << " Go to the East " << std::endl;
while ( ((next=theESNav.east()) != ESDetId(0) && next != strip) ) {
if ( debugLevel_ == pDEBUG ) std::cout << "East: " << n_east << " current strip is " << next << std::endl;
road_2d.push_back(next);
++n_east;
if (n_east == preshSeededNstr_) break;
}
// west road
int n_west= 0;
if ( debugLevel_ == pDEBUG ) std::cout << " Go to the West " << std::endl;
theESNav.home();
while ( ((next=theESNav.west()) != ESDetId(0) && next != strip )) {
if ( debugLevel_ == pDEBUG ) std::cout << "West: " << n_west << " current strip is " << next << std::endl;
road_2d.push_back(next);
++n_west;
if (n_west == preshSeededNstr_) break;
}
if ( debugLevel_ == pDEBUG ) std::cout << "Total number of strips found in the road at 1-st plane is " << n_east+n_west << std::endl;
}
else if (plane == 2) {
// north road
int n_north= 0;
if ( debugLevel_ == pDEBUG ) std::cout << " Go to the North " << std::endl;
while ( ((next=theESNav.north()) != ESDetId(0) && next != strip) ) {
if ( debugLevel_ == pDEBUG ) std::cout << "North: " << n_north << " current strip is " << next << std::endl;
road_2d.push_back(next);
++n_north;
if (n_north == preshSeededNstr_) break;
}
// south road
int n_south= 0;
if ( debugLevel_ == pDEBUG ) std::cout << " Go to the South " << std::endl;
theESNav.home();
while ( ((next=theESNav.south()) != ESDetId(0) && next != strip) ) {
if ( debugLevel_ == pDEBUG ) std::cout << "South: " << n_south << " current strip is " << next << std::endl;
road_2d.push_back(next);
++n_south;
if (n_south == preshSeededNstr_) break;
}
if ( debugLevel_ == pDEBUG ) std::cout << "Total number of strips found in the road at 2-nd plane is " << n_south+n_north << std::endl;
}
else {
if ( debugLevel_ == pDEBUG ) std::cout << " Wrong plane number, null cluster will be returned! " << std::endl;
} // end of if
theESNav.home();
}
| bool PreshowerClusterAlgo::goodStrip | ( | RecHitsMap::iterator | candidate_it | ) |
Definition at line 245 of file PreshowerClusterAlgo.cc.
References gather_cfg::cout, debugLevel_, pDEBUG, preshStripEnergyCut_, rechits_map, and used_s.
Referenced by makeOneCluster().
{
if ( debugLevel_ == pDEBUG ) {
if ( used_s->find(candidate_it->first) != used_s->end())
std::cout << " This strip is in use " << std::endl;
if (candidate_it == rechits_map->end() )
std::cout << " No such a strip in rechits_map " << std::endl;
if (candidate_it->second.energy() <= preshStripEnergyCut_)
std::cout << " Strip energy " << candidate_it->second.energy() <<" is below threshold " << std::endl;
}
// crystal should not be included...
if ( (used_s->find(candidate_it->first) != used_s->end()) || // ...if it already belongs to a cluster
(candidate_it == rechits_map->end() ) || // ...if it corresponds to a hit
(candidate_it->second.energy() <= preshStripEnergyCut_ ) ) // ...if it has a negative or zero energy
{
return false;
}
return true;
}
| reco::PreshowerCluster PreshowerClusterAlgo::makeOneCluster | ( | ESDetId | strip, |
| HitsID * | used_strips, | ||
| RecHitsMap * | rechits_map, | ||
| const CaloSubdetectorGeometry *& | geometry_p, | ||
| CaloSubdetectorTopology *& | topology_p | ||
| ) |
Definition at line 12 of file PreshowerClusterAlgo.cc.
References edm::SortedCollection< T, SORT >::begin(), gather_cfg::cout, CommonMethods::cp(), debugLevel_, CaloNavigator< T >::east(), edm::SortedCollection< T, SORT >::end(), reco::CaloCluster::energy(), reco::CaloCluster::eta(), findRoad(), newFWLiteAna::found, CaloSubdetectorGeometry::getGeometry(), CaloCellGeometry::getPosition(), goodStrip(), CaloNavigator< T >::home(), reco::PreshowerCluster::nhits(), CaloNavigator< T >::north(), pDEBUG, reco::CaloCluster::phi(), pINFO, ESDetId::plane(), pos, reco::CaloCluster::position(), 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().
{
road_2d.clear();
rechits_map = the_rechitsMap_p;
used_s = used_strips;
int plane = strip.plane();
if ( debugLevel_ <= pINFO ) {
std::cout << "Preshower Seeded Algorithm - looking for clusters" << std::endl;
std::cout << "Preshower is intersected at strip " << strip.strip() << ", at plane " << plane << std::endl;
}
// create null-cluster
std::vector< std::pair<DetId,float> > dummy;
Point posi(0,0,0);
if ( debugLevel_ <= pINFO ) std::cout << " Creating a null-cluster" << std::endl;
reco::PreshowerCluster nullcluster=reco::PreshowerCluster(0.,posi,dummy,plane);
if ( strip == ESDetId(0) ) return nullcluster; //works in case of no intersected strip found (e.g. in the Barrel)
// Collection of cluster strips
EcalRecHitCollection clusterRecHits;
// Map of strips for position calculation
RecHitsMap recHits_pos;
//Make a navigator, and set it to the strip cell.
EcalPreshowerNavigator navigator(strip, topology_p);
navigator.setHome(strip);
//search for neighbours in the central road
findRoad(strip,navigator,plane);
if ( debugLevel_ <= pINFO ) std::cout << "Total number of strips in the central road: " << road_2d.size() << std::endl;
if ( plane == 1 ) {
ESDetId strip_north = navigator.north();
findRoad(strip_north,navigator,plane);
navigator.home();
ESDetId strip_south = navigator.south();
findRoad(strip_south,navigator,plane);
navigator.home();
}
if ( plane == 2 ) {
ESDetId strip_east = navigator.east();
findRoad(strip_east,navigator,plane);
navigator.home();
ESDetId strip_west = navigator.west();
findRoad(strip_west,navigator,plane);
navigator.home();
}
if ( debugLevel_ <= pINFO ) std::cout << "Total number of strips in all three roads: " << road_2d.size() << std::endl;
// Start clustering from strip with max Energy in the road
float E_max = 0.;
bool found = false;
RecHitsMap::iterator max_it;
// Loop over strips:
std::vector<ESDetId>::iterator itID;
for (itID = road_2d.begin(); itID != road_2d.end(); itID++) {
if ( debugLevel_ == pDEBUG ) std::cout << " ID = " << *itID << std::endl;
RecHitsMap::iterator strip_it = rechits_map->find(*itID);
//if ( strip_it->second.energy() < 0 ) std::cout << " ##### E = " << strip_it->second.energy() << std::endl;
if(!goodStrip(strip_it)) continue;
if ( debugLevel_ == pDEBUG ) std::cout << " strip is " << ESDetId(strip_it->first) <<" E = " << strip_it->second.energy() <<std::endl;
float E = strip_it->second.energy();
if ( E > E_max) {
E_max = E;
found = true;
max_it = strip_it;
}
}
// no hottest strip found ==> null cluster will be returned
if ( !found ) return nullcluster;
// First, save the hottest strip
clusterRecHits.push_back(max_it->second);
recHits_pos.insert(std::make_pair(max_it->first, max_it->second));
used_s->insert(max_it->first);
if ( debugLevel_ <= pINFO ) {
std::cout << " Central hottest strip " << ESDetId(max_it->first) << " is saved " << std::endl;
std::cout << " with energy E = " << E_max << std::endl;
}
// Find positions of adjacent strips:
ESDetId next, strip_1, strip_2;
navigator.setHome(max_it->first);
ESDetId startES = max_it->first;
if (plane == 1) {
// Save two neighbouring strips to the east
int nadjacents_east = 0;
while ( (next=navigator.east()) != ESDetId(0) && next != startES && nadjacents_east < 2 ) {
++nadjacents_east;
if ( debugLevel_ == pDEBUG ) std::cout << " Adjacent east #" << nadjacents_east <<": "<< next << std::endl;
RecHitsMap::iterator strip_it = rechits_map->find(next);
if(!goodStrip(strip_it)) continue;
// Save strip for clustering if it exists, not already in use, and satisfies an energy threshold
clusterRecHits.push_back(strip_it->second);
// save strip for position calculation
if ( nadjacents_east==1 ) strip_1 = next;
used_s->insert(strip_it->first);
if ( debugLevel_ == pDEBUG ) std::cout << " East adjacent strip # " << nadjacents_east << " is saved with energy E = "
<< strip_it->second.energy() << std::endl;
}
// Save two neighbouring strips to the west
navigator.home();
int nadjacents_west = 0;
while ( (next=navigator.west()) != ESDetId(0) && next != startES && nadjacents_west < 2 ) {
++nadjacents_west;
if ( debugLevel_ == pDEBUG ) std::cout << " Adjacent west #" << nadjacents_west <<": "<< next << std::endl;
RecHitsMap::iterator strip_it = rechits_map->find(next);
if(!goodStrip(strip_it)) continue;
clusterRecHits.push_back(strip_it->second);
if ( nadjacents_west==1 ) strip_2 = next;
used_s->insert(strip_it->first);
if ( debugLevel_ == pDEBUG ) std::cout << " West adjacent strip # " << nadjacents_west << " is saved with energy E = "
<< strip_it->second.energy() << std::endl;
}
}
else if (plane == 2) {
// Save two neighbouring strips to the north
int nadjacents_north = 0;
while ( (next=navigator.north()) != ESDetId(0) && next != startES && nadjacents_north < 2 ) {
++nadjacents_north;
if ( debugLevel_ == pDEBUG ) std::cout << " Adjacent north #" << nadjacents_north <<": "<< next << std::endl;
RecHitsMap::iterator strip_it = rechits_map->find(next);
if(!goodStrip(strip_it)) continue;
clusterRecHits.push_back(strip_it->second);
if ( nadjacents_north==1 ) strip_1 = next;
used_s->insert(strip_it->first);
if ( debugLevel_ == pDEBUG ) std::cout << " North adjacent strip # " << nadjacents_north << " is saved with energy E = "
<< strip_it->second.energy() << std::endl;
}
// Save two neighbouring strips to the south
navigator.home();
int nadjacents_south = 0;
while ( (next=navigator.south()) != ESDetId(0) && next != startES && nadjacents_south < 2 ) {
++nadjacents_south;
if ( debugLevel_ == pDEBUG ) std::cout << " Adjacent south #" << nadjacents_south <<": "<< next << std::endl;
RecHitsMap::iterator strip_it = rechits_map->find(next);
if(!goodStrip(strip_it)) continue;
clusterRecHits.push_back(strip_it->second);
if ( nadjacents_south==1 ) strip_2 = next;
used_s->insert(strip_it->first);
if ( debugLevel_ == pDEBUG ) std::cout << " South adjacent strip # " << nadjacents_south << " is saved with energy E = "
<< strip_it->second.energy() << std::endl;
}
}
else {
std::cout << " Wrong plane number" << plane <<", null cluster will be returned! " << std::endl;
return nullcluster;
} // end of if
if ( debugLevel_ <=pINFO ) std::cout << " Total size of clusterRecHits is " << clusterRecHits.size() << std::endl;
if ( debugLevel_ <=pINFO ) std::cout << " Two adjacent strips for position calculation are: "
<< strip_1 <<" and " << strip_2 << std::endl;
// strips for position calculation
RecHitsMap::iterator strip_it1, strip_it2;
if ( strip_1 != ESDetId(0)) {
strip_it1 = rechits_map->find(strip_1);
recHits_pos.insert(std::make_pair(strip_it1->first, strip_it1->second));
}
if ( strip_2 != ESDetId(0) ) {
strip_it2 = rechits_map->find(strip_2);
recHits_pos.insert(std::make_pair(strip_it2->first, strip_it2->second));
}
RecHitsMap::iterator cp;
double energy_pos = 0;
double x_pos = 0;
double y_pos = 0;
double z_pos = 0;
for (cp = recHits_pos.begin(); cp!=recHits_pos.end(); cp++ ) {
double E = cp->second.energy();
energy_pos += E;
const CaloCellGeometry *thisCell = geometry_p->getGeometry(cp->first);
GlobalPoint position = thisCell->getPosition();
x_pos += E * position.x();
y_pos += E * position.y();
z_pos += E * position.z();
}
if(energy_pos>0.) {
x_pos /= energy_pos;
y_pos /= energy_pos;
z_pos /= energy_pos;
}
Point pos(x_pos,y_pos,z_pos);
if ( debugLevel_ == pDEBUG ) std::cout << " ES Cluster position = " << "(" << x_pos <<","<< y_pos <<","<< z_pos <<")"<< std::endl;
EcalRecHitCollection::iterator it;
double Eclust = 0;
if ( debugLevel_ == pINFO ) std::cout << "The found ES cluster strips: " << std::endl;
std::vector<std::pair<DetId,float > > usedHits;
for (it=clusterRecHits.begin(); it != clusterRecHits.end(); it++) {
Eclust += it->energy();
usedHits.push_back(std::pair<DetId,float > (it->id(),1.));
if ( debugLevel_ == pINFO ) std::cout << ESDetId(it->id()) <<", E = " << it->energy()<<"; ";
}
if ( debugLevel_ == pINFO ) std::cout << std::endl;
// ES cluster is created from vector clusterRecHits
reco::PreshowerCluster cluster=reco::PreshowerCluster(Eclust,pos,usedHits,plane);
if ( debugLevel_ <= pINFO ) {
std::cout << " ES Cluster is created with " << std::endl;
std::cout << " energy = " << cluster.energy() << std::endl;
std::cout << " (eta,phi) = " << "("<<cluster.eta()<<", "<<cluster.phi()<<")"<< std::endl;
std::cout << " nhits = " << cluster.nhits() << std::endl;
std::cout << " radius = " << cluster.position().r() << std::endl;
std::cout << " (x,y,z) = " << "(" << cluster.x() <<", "<< cluster.y() <<", "<< cluster.z()<<")"<< std::endl;
}
used_strips = used_s;
// return the cluster if its energy is greater a threshold
if( cluster.energy() > preshClusterEnergyCut_ )
return cluster;
else
return nullcluster;
}
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().
1.7.3