File indexing completed on 2023-03-17 11:17:19
0001 #include "RecoEcal/EgammaClusterAlgos/interface/IslandClusterAlgo.h"
0002
0003
0004 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0005 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0006 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0007 #include "Geometry/CaloTopology/interface/EcalEndcapTopology.h"
0008 #include "Geometry/CaloTopology/interface/EcalBarrelTopology.h"
0009
0010 #include "RecoEcal/EgammaCoreTools/interface/PositionCalc.h"
0011 #include "RecoEcal/EgammaCoreTools/interface/ClusterEtLess.h"
0012
0013
0014 #include "DataFormats/CaloRecHit/interface/CaloID.h"
0015
0016
0017 std::vector<reco::BasicCluster> IslandClusterAlgo::makeClusters(const EcalRecHitCollection *hits,
0018 const CaloSubdetectorGeometry *geometry_p,
0019 const CaloSubdetectorTopology *topology_p,
0020 const CaloSubdetectorGeometry *geometryES_p,
0021 EcalPart ecalPart,
0022 bool regional,
0023 const std::vector<RectangularEtaPhiRegion> ®ions) {
0024 seeds.clear();
0025 used_s.clear();
0026 clusters_v.clear();
0027
0028 recHits_ = hits;
0029
0030 double threshold = 0;
0031 std::string ecalPart_string;
0032 if (ecalPart == endcap) {
0033 threshold = ecalEndcapSeedThreshold;
0034 ecalPart_string = "EndCap";
0035 v_chstatusSeed_ = v_chstatusSeed_Endcap_;
0036 v_chstatus_ = v_chstatus_Endcap_;
0037 }
0038 if (ecalPart == barrel) {
0039 threshold = ecalBarrelSeedThreshold;
0040 ecalPart_string = "Barrel";
0041 v_chstatusSeed_ = v_chstatusSeed_Barrel_;
0042 v_chstatus_ = v_chstatus_Barrel_;
0043 }
0044
0045 if (verbosity < pINFO) {
0046 std::cout << "-------------------------------------------------------------" << std::endl;
0047 std::cout << "Island algorithm invoked for ECAL" << ecalPart_string << std::endl;
0048 std::cout << "Looking for seeds, energy threshold used = " << threshold << " GeV" << std::endl;
0049 }
0050
0051 int nregions = 0;
0052 if (regional)
0053 nregions = regions.size();
0054
0055 if (!regional || nregions) {
0056 EcalRecHitCollection::const_iterator it;
0057 for (it = hits->begin(); it != hits->end(); it++) {
0058 double energy = it->energy();
0059 if (energy < threshold)
0060 continue;
0061
0062
0063 if (!it->checkFlag(EcalRecHit::kGood)) {
0064 if (it->checkFlags(v_chstatus_) || it->checkFlags(v_chstatusSeed_)) {
0065 continue;
0066 }
0067 }
0068
0069 auto thisCell = geometry_p->getGeometry(it->id());
0070 auto const &position = thisCell->getPosition();
0071
0072
0073
0074 bool withinRegion = false;
0075 if (regional) {
0076 std::vector<RectangularEtaPhiRegion>::const_iterator region;
0077 for (region = regions.begin(); region != regions.end(); region++) {
0078 if (region->inRegion(thisCell->etaPos(), thisCell->phiPos())) {
0079 withinRegion = true;
0080 break;
0081 }
0082 }
0083 }
0084
0085 if (!regional || withinRegion) {
0086 float ET = it->energy() * position.basicVector().unit().perp();
0087 if (ET > threshold)
0088 seeds.push_back(*it);
0089 }
0090 }
0091 }
0092
0093 sort(seeds.begin(), seeds.end(), [](auto const &x, auto const &y) { return x.energy() > y.energy(); });
0094
0095 if (verbosity < pINFO) {
0096 std::cout << "Total number of seeds found in event = " << seeds.size() << std::endl;
0097 }
0098
0099 mainSearch(hits, geometry_p, topology_p, geometryES_p, ecalPart);
0100 sort(clusters_v.rbegin(), clusters_v.rend(), isClusterEtLess);
0101
0102 if (verbosity < pINFO) {
0103 std::cout << "---------- end of main search. clusters have been sorted ----" << std::endl;
0104 }
0105
0106 return clusters_v;
0107 }
0108
0109 void IslandClusterAlgo::mainSearch(const EcalRecHitCollection *hits,
0110 const CaloSubdetectorGeometry *geometry_p,
0111 const CaloSubdetectorTopology *topology_p,
0112 const CaloSubdetectorGeometry *geometryES_p,
0113 EcalPart ecalPart) {
0114 if (verbosity < pINFO) {
0115 std::cout << "Building clusters............" << std::endl;
0116 }
0117
0118
0119 std::vector<EcalRecHit>::iterator it;
0120 for (it = seeds.begin(); it != seeds.end(); it++) {
0121
0122 if (used_s.find(it->id()) != used_s.end()) {
0123 if (it == seeds.begin()) {
0124 if (verbosity < pINFO) {
0125 std::cout << "##############################################################" << std::endl;
0126 std::cout << "DEBUG ALERT: Highest energy seed already belongs to a cluster!" << std::endl;
0127 std::cout << "##############################################################" << std::endl;
0128 }
0129 }
0130 continue;
0131 }
0132
0133
0134 current_v.clear();
0135
0136 current_v.push_back(std::pair<DetId, float>(it->id(), 1.));
0137 used_s.insert(it->id());
0138
0139
0140 CaloNavigator<DetId> navigator(it->id(), topology_p);
0141
0142 searchNorth(navigator);
0143 navigator.home();
0144 searchSouth(navigator);
0145 navigator.home();
0146 searchWest(navigator, topology_p);
0147 navigator.home();
0148 searchEast(navigator, topology_p);
0149
0150 makeCluster(hits, geometry_p, geometryES_p);
0151 }
0152 }
0153
0154 void IslandClusterAlgo::searchNorth(const CaloNavigator<DetId> &navigator) {
0155 DetId southern = navigator.pos();
0156
0157 DetId northern = navigator.north();
0158 if (northern == DetId(0))
0159 return;
0160
0161 if (used_s.find(northern) != used_s.end())
0162 return;
0163
0164 EcalRecHitCollection::const_iterator southern_it = recHits_->find(southern);
0165 EcalRecHitCollection::const_iterator northern_it = recHits_->find(northern);
0166
0167 if (shouldBeAdded(northern_it, southern_it)) {
0168 current_v.push_back(std::pair<DetId, float>(northern, 1.));
0169 used_s.insert(northern);
0170 searchNorth(navigator);
0171 }
0172 }
0173
0174 void IslandClusterAlgo::searchSouth(const CaloNavigator<DetId> &navigator) {
0175 DetId northern = navigator.pos();
0176
0177 DetId southern = navigator.south();
0178 if (southern == DetId(0))
0179 return;
0180 if (used_s.find(southern) != used_s.end())
0181 return;
0182
0183 EcalRecHitCollection::const_iterator northern_it = recHits_->find(northern);
0184 EcalRecHitCollection::const_iterator southern_it = recHits_->find(southern);
0185
0186 if (shouldBeAdded(southern_it, northern_it)) {
0187 current_v.push_back(std::pair<DetId, float>(southern, 1.));
0188 used_s.insert(southern);
0189 searchSouth(navigator);
0190 }
0191 }
0192
0193 void IslandClusterAlgo::searchWest(const CaloNavigator<DetId> &navigator, const CaloSubdetectorTopology *topology) {
0194 DetId eastern = navigator.pos();
0195 EcalRecHitCollection::const_iterator eastern_it = recHits_->find(eastern);
0196
0197 DetId western = navigator.west();
0198 if (western == DetId(0))
0199 return;
0200 EcalRecHitCollection::const_iterator western_it = recHits_->find(western);
0201
0202 if (shouldBeAdded(western_it, eastern_it)) {
0203 CaloNavigator<DetId> nsNavigator(western, topology);
0204
0205 searchNorth(nsNavigator);
0206 nsNavigator.home();
0207 searchSouth(nsNavigator);
0208 nsNavigator.home();
0209 searchWest(navigator, topology);
0210
0211 current_v.push_back(std::pair<DetId, float>(western, 1.));
0212 used_s.insert(western);
0213 }
0214 }
0215
0216 void IslandClusterAlgo::searchEast(const CaloNavigator<DetId> &navigator, const CaloSubdetectorTopology *topology) {
0217 DetId western = navigator.pos();
0218 EcalRecHitCollection::const_iterator western_it = recHits_->find(western);
0219
0220 DetId eastern = navigator.east();
0221 if (eastern == DetId(0))
0222 return;
0223 EcalRecHitCollection::const_iterator eastern_it = recHits_->find(eastern);
0224
0225 if (shouldBeAdded(eastern_it, western_it)) {
0226 CaloNavigator<DetId> nsNavigator(eastern, topology);
0227
0228 searchNorth(nsNavigator);
0229 nsNavigator.home();
0230 searchSouth(nsNavigator);
0231 nsNavigator.home();
0232 searchEast(navigator, topology);
0233
0234 current_v.push_back(std::pair<DetId, float>(eastern, 1.));
0235 used_s.insert(eastern);
0236 }
0237 }
0238
0239
0240 bool IslandClusterAlgo::shouldBeAdded(EcalRecHitCollection::const_iterator candidate_it,
0241 EcalRecHitCollection::const_iterator previous_it) {
0242
0243 if ((candidate_it == recHits_->end()) ||
0244 (used_s.find(candidate_it->id()) != used_s.end()) ||
0245 (candidate_it->energy() <= 0) ||
0246 (candidate_it->energy() > previous_it->energy()) ||
0247 (!(candidate_it->checkFlag(EcalRecHit::kGood)) && candidate_it->checkFlags(v_chstatus_))) {
0248 return false;
0249 }
0250 return true;
0251 }
0252
0253 void IslandClusterAlgo::makeCluster(const EcalRecHitCollection *hits,
0254 const CaloSubdetectorGeometry *geometry,
0255 const CaloSubdetectorGeometry *geometryES) {
0256 double energy = 0;
0257 reco::CaloID caloID;
0258
0259 Point position;
0260 position = posCalculator_.Calculate_Location(current_v, hits, geometry, geometryES);
0261
0262 std::vector<std::pair<DetId, float> >::iterator it;
0263 for (it = current_v.begin(); it != current_v.end(); it++) {
0264 EcalRecHitCollection::const_iterator itt = hits->find((*it).first);
0265 EcalRecHit hit_p = *itt;
0266 if ((*it).first.subdetId() == EcalBarrel) {
0267 caloID = reco::CaloID::DET_ECAL_BARREL;
0268 } else {
0269 caloID = reco::CaloID::DET_ECAL_ENDCAP;
0270 }
0271
0272
0273 energy += hit_p.energy();
0274
0275
0276
0277
0278
0279 }
0280
0281 if (verbosity < pINFO) {
0282 std::cout << "******** NEW CLUSTER ********" << std::endl;
0283 std::cout << "No. of crystals = " << current_v.size() << std::endl;
0284 std::cout << " Energy = " << energy << std::endl;
0285 std::cout << " Phi = " << position.phi() << std::endl;
0286 std::cout << " Eta = " << position.eta() << std::endl;
0287 std::cout << "*****************************" << std::endl;
0288 }
0289 clusters_v.push_back(reco::BasicCluster(energy, position, caloID, current_v, reco::CaloCluster::island));
0290 }