|
|
|
#include <PreshowerClusterAlgo.h>
Public Types | |
| 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) | |
| PreshowerClusterAlgo () | |
| ~PreshowerClusterAlgo () | |
Private Attributes | |
| 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 33 of file PreshowerClusterAlgo.h.
Definition at line 30 of file PreshowerClusterAlgo.h.
| typedef std::map<DetId, EcalRecHit> PreshowerClusterAlgo::RecHitsMap |
Definition at line 32 of file PreshowerClusterAlgo.h.
| PreshowerClusterAlgo::PreshowerClusterAlgo | ( | ) | [inline] |
Definition at line 35 of file PreshowerClusterAlgo.h.
: preshStripEnergyCut_(0.), preshClusterEnergyCut_(0.), preshSeededNstr_(15) {}
| PreshowerClusterAlgo::PreshowerClusterAlgo | ( | double | stripEnergyCut, |
| double | clusterEnergyCut, | ||
| int | nStripCut | ||
| ) | [inline] |
Definition at line 39 of file PreshowerClusterAlgo.h.
: preshStripEnergyCut_(stripEnergyCut), preshClusterEnergyCut_(clusterEnergyCut), preshSeededNstr_(nStripCut) {}
| PreshowerClusterAlgo::~PreshowerClusterAlgo | ( | ) | [inline] |
Definition at line 43 of file PreshowerClusterAlgo.h.
{};
| void PreshowerClusterAlgo::findRoad | ( | ESDetId | strip, |
| EcalPreshowerNavigator | theESNav, | ||
| int | plane | ||
| ) |
Definition at line 268 of file PreshowerClusterAlgo.cc.
References CaloNavigator< T >::east(), CaloNavigator< T >::home(), LogTrace, GetRecoTauVFromDQM_MC_cff::next, CaloNavigator< T >::north(), 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);
LogTrace("PreShowerClusterAlgo") << "findRoad starts from strip" << strip;
if (plane == 1) {
// east road
int n_east= 0;
LogTrace("PreShowerClusterAlgo") << " Go to the East ";
while ( ((next=theESNav.east()) != ESDetId(0) && next != strip) ) {
LogTrace("PreShowerClusterAlgo") << "East:" << n_east << "current strip is"<< next;
road_2d.push_back(next);
++n_east;
if (n_east == preshSeededNstr_) break;
}
// west road
int n_west= 0;
LogTrace("PreShowerClusterAlgo") << "Go to the West" ;
theESNav.home();
while ( ((next=theESNav.west()) != ESDetId(0) && next != strip )) {
LogTrace("PreShowerClusterAlgo") << "West: " << n_west << "current strip is" << next ;
road_2d.push_back(next);
++n_west;
if (n_west == preshSeededNstr_) break;
}
LogTrace("PreShowerClusterAlgo") << "Total number of strips found in the road at 1-st plane is" << n_east+n_west ;
}
else if (plane == 2) {
// north road
int n_north= 0;
LogTrace("PreShowerClusterAlgo") << "Go to the North";
while ( ((next=theESNav.north()) != ESDetId(0) && next != strip) ) {
LogTrace("PreShowerClusterAlgo") << "North:" << n_north << "current strip is" << next ;
road_2d.push_back(next);
++n_north;
if (n_north == preshSeededNstr_) break;
}
// south road
int n_south= 0;
LogTrace("PreShowerClusterAlgo") << "Go to the South";
theESNav.home();
while ( ((next=theESNav.south()) != ESDetId(0) && next != strip) ) {
LogTrace("PreShowerClusterAlgo") << "South:" << n_south << "current strip is" << next ;
road_2d.push_back(next);
++n_south;
if (n_south == preshSeededNstr_) break;
}
LogTrace("PreShowerClusterAlgo") << "Total number of strips found in the road at 2-nd plane is" << n_south+n_north;
}
else {
LogTrace("PreShowerClusterAlgo") << " Wrong plane number, null cluster will be returned!";
} // end of if
theESNav.home();
}
| bool PreshowerClusterAlgo::goodStrip | ( | RecHitsMap::iterator | candidate_it | ) |
Definition at line 246 of file PreshowerClusterAlgo.cc.
References LogTrace, preshStripEnergyCut_, rechits_map, and used_s.
Referenced by makeOneCluster().
{
if ( used_s->find(candidate_it->first) != used_s->end())
LogTrace("PreShowerClusterAlgo") << " This strip is in use";
if (candidate_it == rechits_map->end() )
LogTrace("PreShowerClusterAlgo") << " No such a strip in rechits_map";
if (candidate_it->second.energy() <= preshStripEnergyCut_)
LogTrace("PreShowerClusterAlgo") << "Strip energy" << candidate_it->second.energy() <<"is below threshold";
// 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 13 of file PreshowerClusterAlgo.cc.
References edm::SortedCollection< T, SORT >::begin(), CommonMethods::cp(), 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(), LogTrace, GetRecoTauVFromDQM_MC_cff::next, reco::PreshowerCluster::nhits(), CaloNavigator< T >::north(), reco::CaloCluster::phi(), 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();
LogTrace("PreShowerClusterAlgo") << "Preshower Seeded Algorithm - looking for clusters";
LogTrace("PreShowerClusterAlgo")<< "Preshower is intersected at strip" << strip.strip() << ",at plane" << plane ;
// create null-cluster
std::vector< std::pair<DetId,float> > dummy;
Point posi(0,0,0);
LogTrace("PreShowerClusterAlgo") << " Creating a null-cluster" ;
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);
LogTrace("PreShowerClusterAlgo") << "Total number of strips in the central road:" << road_2d.size() ;
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();
}
LogTrace("PreShowerClusterAlgo") << "Total number of strips in all three roads:" << road_2d.size() ;
// 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++) {
LogTrace("PreShowerClusterAlgo") << "ID ="<<*itID ;
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;
LogTrace("PreShowerClusterAlgo") << " strip is " << ESDetId(strip_it->first) <<"E ="<< strip_it->second.energy();
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);
LogTrace("PreShowerClusterAlgo") << "Central hottest strip" << ESDetId(max_it->first) << "is saved with energy E =" << E_max ;
// 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;
LogTrace("PreShowerClusterAlgo") << "Adjacent east #" << nadjacents_east <<":"<< next ;
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);
LogTrace("PreShowerClusterAlgo") << "East adjacent strip #" << nadjacents_east << "is saved with energy E =" << strip_it->second.energy() ;
}
// 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;
LogTrace("PreShowerClusterAlgo") << "Adjacent west #" << nadjacents_west <<":"<< next ;
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);
LogTrace("PreShowerClusterAlgo") << "West adjacent strip #" << nadjacents_west << "is saved with energy E =" << strip_it->second.energy();
}
}
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;
LogTrace("PreShowerClusterAlgo") << "Adjacent north #" << nadjacents_north <<":"<< next ;
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);
LogTrace("PreShowerClusterAlgo") << "North adjacent strip #" << nadjacents_north << "is saved with energy E =" << strip_it->second.energy() ;
}
// 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;
LogTrace("PreShowerClusterAlgo") << "Adjacent south #" << nadjacents_south <<":"<< next ;
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);
LogTrace("PreShowerClusterAlgo") << "South adjacent strip #" << nadjacents_south << "is saved with energy E =" << strip_it->second.energy() ;
}
}
else {
LogTrace("PreShowerClusterAlgo") << " Wrong plane number" << plane <<", null cluster will be returned! " << std::endl;
return nullcluster;
} // end of if
LogTrace("PreShowerClusterAlgo") << "Total size of clusterRecHits is" << clusterRecHits.size();
LogTrace("PreShowerClusterAlgo") << "Two adjacent strips for position calculation are:" << strip_1 <<"and" << strip_2;
// 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);
LogTrace("PreShowerClusterAlgo") << "ES Cluster position =" << "("<< x_pos <<","<< y_pos <<","<< z_pos <<")";
EcalRecHitCollection::iterator it;
double Eclust = 0;
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.));
}
// ES cluster is created from vector clusterRecHits
reco::PreshowerCluster cluster=reco::PreshowerCluster(Eclust,pos,usedHits,plane);
LogTrace("PreShowerClusterAlgo") << " ES Cluster is created with:" ;
LogTrace("PreShowerClusterAlgo") << " energy =" << cluster.energy();
LogTrace("PreShowerClusterAlgo") << " (eta,phi) =" << "("<<cluster.eta()<<","<<cluster.phi()<<")";
LogTrace("PreShowerClusterAlgo") << " nhits =" << cluster.nhits();
LogTrace("PreShowerClusterAlgo") << " radius =" << cluster.position().r();
LogTrace("PreShowerClusterAlgo") << " (x,y,z) =" << "(" << cluster.x() <<", "<< cluster.y() <<","<< cluster.z()<<")" ;
used_strips = used_s;
// return the cluster if its energy is greater a threshold
if( cluster.energy() > preshClusterEnergyCut_ )
return cluster;
else
return nullcluster;
}
double PreshowerClusterAlgo::preshClusterEnergyCut_ [private] |
Definition at line 58 of file PreshowerClusterAlgo.h.
Referenced by makeOneCluster().
int PreshowerClusterAlgo::preshSeededNstr_ [private] |
Definition at line 59 of file PreshowerClusterAlgo.h.
Referenced by findRoad().
double PreshowerClusterAlgo::preshStripEnergyCut_ [private] |
Definition at line 57 of file PreshowerClusterAlgo.h.
Referenced by goodStrip().
RecHitsMap* PreshowerClusterAlgo::rechits_map [private] |
Definition at line 65 of file PreshowerClusterAlgo.h.
Referenced by goodStrip(), and makeOneCluster().
std::vector<ESDetId> PreshowerClusterAlgo::road_2d [private] |
Definition at line 62 of file PreshowerClusterAlgo.h.
Referenced by findRoad(), and makeOneCluster().
HitsID* PreshowerClusterAlgo::used_s [private] |
Definition at line 68 of file PreshowerClusterAlgo.h.
Referenced by goodStrip(), and makeOneCluster().
1.7.3