16 template <
class T1,
class T2,
typename Comp=std::less<T1> >
struct PairSortByFirst {
18 bool operator()(
const std::pair<T1,T2>& lhs,
const std::pair<T1,T2>&rhs){
return comp(lhs.first,rhs.first);}
19 bool operator()(
const T1& lhs,
const std::pair<T1,T2>&rhs){
return comp(lhs,rhs.first);}
20 bool operator()(
const std::pair<T1,T2>& lhs,
const T1 &rhs){
return comp(lhs.first,rhs);}
26 if(lhs.
ieta()<rhs.
ieta())
return true;
27 else if(lhs.
ieta()>rhs.
ieta())
return false;
37 if(lhs.
ix()<rhs.
ix())
return true;
38 else if(lhs.
ix()>rhs.
ix())
return false;
39 else return lhs.
iy()<rhs.
iy();
53 const std::vector<EcalEtaPhiRegion>& regions
69 threshold = ecalEndcapSeedThreshold;
70 ecalPart_string =
"EndCap";
75 threshold = ecalBarrelSeedThreshold;
76 ecalPart_string =
"Barrel";
79 LogTrace(
"EcalClusters") <<
"-------------------------------------------------------------";
80 LogTrace(
"EcalClusters") <<
"Island algorithm invoked for ECAL" << ecalPart_string ;
81 LogTrace(
"EcalClusters") <<
"Looking for seeds, energy threshold used = " << threshold <<
" GeV";
85 if(regional) nregions=regions.size();
87 if(!regional || nregions) {
90 for(it = hits->
begin(); it != hits->
end(); it++)
92 double energy = it->energy();
93 if (energy < threshold)
continue;
98 bool withinRegion =
false;
100 std::vector<EcalEtaPhiRegion>::const_iterator region;
101 for (region=regions.begin(); region!=regions.end(); region++) {
102 if (region->inRegion(thisCell->etaPos(),thisCell->phiPos())) {
109 if (!regional || withinRegion) {
110 float ET = it->energy() * thisCell->getPosition().basicVector().unit().perp();
111 if (ET > threshold) seeds.push_back(*it);
119 LogTrace(
"EcalClusters") <<
"Total number of seeds found in event = " << seeds.size();
122 mainSearch(hits, geometry_p, topology_p, geometryES_p);
123 sort(clusters_v.rbegin(), clusters_v.rend(),
ClusterEtLess());
125 LogTrace(
"EcalClusters") <<
"---------- end of main search. clusters have been sorted ----";
141 LogTrace(
"EcalClusters") <<
"Building clusters............";
144 std::vector<EcalRecHit>::iterator it;
145 for (it = seeds.begin(); it != seeds.end(); it++)
150 bool usedButCanSeed =
false;
151 if (canSeed_s.find(it->id()) != canSeed_s.end()) usedButCanSeed =
true;
155 if (it->checkFlags( v_chstatus_ ))
continue;
159 if ((used_s.find(it->id()) != used_s.end()) && (usedButCanSeed ==
false))
161 if (it == seeds.begin())
163 LogTrace(
"EcalClusters") <<
"##############################################################" ;
164 LogTrace(
"EcalClusters") <<
"DEBUG ALERT: Highest energy seed already belongs to a cluster!";
165 LogTrace(
"EcalClusters") <<
"##############################################################";
185 bool localMaxima = checkMaxima(
navigator, hits);
191 prepareCluster(
navigator, hits, geometry_p);
196 if (!current_v.empty())
198 makeCluster(hits, geometry_p, geometryES_p, seedIt, usedButCanSeed);
204 if(reassignSeedCrysToClusterItSeeds_){
208 for(
size_t clusNr=0;clusNr<protoClusters_.size();clusNr++){
209 if(!protoClusters_[clusNr].containsSeed()){
210 const EcalRecHit& seedHit =protoClusters_[clusNr].seed();
211 typedef std::vector<std::pair<DetId,int> >::iterator It;
214 if(result.first!=result.second) protoClusters_[result.first->second].removeHit(seedHit);
215 protoClusters_[clusNr].addSeed();
223 for(
size_t clusNr=0;clusNr<protoClusters_.size();clusNr++){
226 position = posCalculator_.Calculate_Location(protoCluster.
hits(),
hits, geometry_p, geometryES_p);
231 protoClusters_.clear();
232 whichClusCrysBelongsTo_.clear();
246 position = posCalculator_.Calculate_Location(current_v, hits, geometry, geometryES);
248 std::vector<std::pair<DetId, float> >::iterator it;
249 for (it = current_v.begin(); it != current_v.end(); it++)
264 LogTrace(
"EcalClusters") <<
"******** NEW CLUSTER ********";
265 LogTrace(
"EcalClusters") <<
"No. of crystals = " << current_v.size();
266 LogTrace(
"EcalClusters") <<
" Energy = " << energy ;
267 LogTrace(
"EcalClusters") <<
" Phi = " << position.phi();
268 LogTrace(
"EcalClusters") <<
" Eta " << position.eta();
269 LogTrace(
"EcalClusters") <<
"*****************************";
274 double seedEnergy = seedIt->energy();
275 if ((seedOutside && energy>=0) || (!seedOutside && energy >= seedEnergy))
277 if(reassignSeedCrysToClusterItSeeds_){
278 for(
size_t hitNr=0;hitNr<current_v.size();hitNr++) whichClusCrysBelongsTo_.push_back(std::pair<DetId,int>(current_v[hitNr].first,protoClusters_.size()));
287 std::vector<std::pair<DetId, float> >::iterator iter;
288 for (iter = current_v.begin(); iter != current_v.end(); iter++)
290 used_s.erase(iter->first);
303 double seedEnergy = seedHit->energy();
305 std::vector<DetId> swissCrossVec;
306 swissCrossVec.clear();
308 swissCrossVec.push_back(navigator.
west());
310 swissCrossVec.push_back(navigator.
east());
312 swissCrossVec.push_back(navigator.
north());
314 swissCrossVec.push_back(navigator.
south());
317 std::vector<DetId>::const_iterator detItr;
318 for (
unsigned int i = 0;
i < swissCrossVec.size(); ++
i)
322 thisHit = recHits_->find(swissCrossVec[
i]);
325 if ((swissCrossVec[i] ==
DetId(0)) || thisHit == recHits_->end())
continue;
329 uint32_t rhFlag = thisHit->recoFlag();
330 std::vector<int>::const_iterator vit =
std::find(v_chstatus_.begin(), v_chstatus_.end(), rhFlag);
331 if (vit != v_chstatus_.end())
continue;
335 if (thisHit->energy() > seedEnergy)
352 std::set<DetId>::iterator setItr;
359 for (
int dx = -2;
dx < 3; ++
dx)
361 for (
int dy = -2;
dy < 3; ++
dy)
380 canSeed_s.insert(thisDet);
386 setItr = canSeed_s.find(thisDet);
387 if (setItr != canSeed_s.end())
390 canSeed_s.erase(setItr);
409 if ((thisIt != recHits_->end()) && (thisIt->id() !=
DetId(0)))
411 if ((used_s.find(thisIt->id()) == used_s.end()))
414 current_v.push_back( std::pair<DetId, float>(det, 1.) );
423 std::vector<std::pair<DetId,float> >::iterator hitPos;
424 for(hitPos=hits_.begin();hitPos<hits_.end();hitPos++){
425 if(hitToRM.
id()==hitPos->first)
break;
427 if(hitPos!=hits_.end()){
429 energy_-=hitToRM.
energy();
438 typedef std::vector<std::pair<DetId,float> >::iterator It;
439 std::pair<It,It> hitPos;
447 if(hitPos.first==hitPos.second){
448 hits_.insert(hitPos.first,std::pair<DetId,float>(seed_.id(),1.));
449 energy_+=seed_.energy();
459 for(
size_t hitNr=0;hitNr<hits_.size();hitNr++){
460 if(seed_.id()==hits_[hitNr].first)
return true;
void makeCluster(const EcalRecHitCollection *hits, const CaloSubdetectorGeometry *geometry_p, const CaloSubdetectorGeometry *geometryES_p, const EcalRecHitCollection::const_iterator &seedIt, bool seedOutside)
std::vector< reco::BasicCluster > makeClusters(const EcalRecHitCollection *hits, const CaloSubdetectorGeometry *geometry, const CaloSubdetectorTopology *topology_p, const CaloSubdetectorGeometry *geometryES_p, reco::CaloID::Detectors detector, bool regional=false, const std::vector< EcalEtaPhiRegion > ®ions=std::vector< EcalEtaPhiRegion >())
std::vector< EcalRecHit >::const_iterator const_iterator
const std::vector< std::pair< DetId, float > > & hits() const
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
int iphi() const
get the crystal iphi
void mainSearch(const EcalRecHitCollection *hits, const CaloSubdetectorGeometry *geometry_p, const CaloSubdetectorTopology *topology_p, const CaloSubdetectorGeometry *geometryES_p)
T offsetBy(int deltaX, int deltaY) const
Free movement of arbitray steps.
void prepareCluster(CaloNavigator< DetId > &navigator, const EcalRecHitCollection *hits, const CaloSubdetectorGeometry *geometry)
bool operator()(const std::pair< T1, T2 > &lhs, const T1 &rhs)
T west() const
move the navigator west
void addCrystal(const DetId &det)
Abs< T >::type abs(const T &t)
bool operator()(const EEDetId &lhs, const EEDetId &rhs)
T south() const
move the navigator south
int ieta() const
get the crystal ieta
T pos() const
get the current position
bool operator()(const T1 &lhs, const T1 &rhs)
const_iterator end() const
bool isSeedCrysInHits_() const
bool removeHit(const EcalRecHit &hitToRM)
T east() const
move the navigator east
void home() const
move the navigator back to the starting point
DetId id() const
get the id
bool operator()(const EBDetId &lhs, const EBDetId &rhs)
bool operator()(const std::pair< T1, T2 > &lhs, const std::pair< T1, T2 > &rhs)
virtual std::shared_ptr< const CaloCellGeometry > getGeometry(const DetId &id) const
Get the cell geometry of a given detector id. Should return false if not found.
bool checkMaxima(CaloNavigator< DetId > &navigator, const EcalRecHitCollection *hits)
bool operator()(const T1 &lhs, const std::pair< T1, T2 > &rhs)
math::XYZPoint Point
point in the space
iterator find(key_type k)
static int position[264][3]
T north() const
move the navigator north
const EcalRecHit & seed() const
const_iterator begin() const