EgammaClusterAlgos/src/IslandClusterAlgo.cc

0001 #include "RecoEcal/EgammaClusterAlgos/interface/IslandClusterAlgo.h"
0002
0003 // Geometry
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 // Return a vector of clusters from a collection of EcalRecHits:
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> &regions) {
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;  // need to check to see if this line is useful!
0061
0062       // avoid seeding for anomalous channels
0063       if (!it->checkFlag(EcalRecHit::kGood)) {  // if rechit is good, no need for further checks
0064         if (it->checkFlags(v_chstatus_) || it->checkFlags(v_chstatusSeed_)) {
0065           continue;  // the recHit has to be excluded from seeding
0066         }
0067       }
0068
0069       auto thisCell = geometry_p->getGeometry(it->id());
0070       auto const &position = thisCell->getPosition();
0071
0072       // Require that RecHit is within clustering region in case
0073       // of regional reconstruction
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   // Loop over seeds:
0119   std::vector<EcalRecHit>::iterator it;
0120   for (it = seeds.begin(); it != seeds.end(); it++) {
0121     // make sure the current seed does not belong to a cluster already.
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     // clear the vector of hits in current cluster
0134     current_v.clear();
0135
0136     current_v.push_back(std::pair<DetId, float>(it->id(), 1.));  // by default hit energy fractions are set at 1.
0137     used_s.insert(it->id());
0138
0139     // Create a navigator at the seed
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;  // This means that we went off the ECAL!
0160   // if the crystal to the north belongs to another cluster return
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.));  // by default hit energy fractions are set at 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;  // This means that we went off the ECAL!
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.));  // by default hit energy fractions are set at 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;  // This means that we went off the ECAL!
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.));  // by default hit energy fractions are set at 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;  // This means that we went off the ECAL!
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.));  // by default hit energy fractions are set at 1.
0235     used_s.insert(eastern);
0236   }
0237 }
0238
0239 // returns true if the candidate crystal fulfills the requirements to be added to the cluster:
0240 bool IslandClusterAlgo::shouldBeAdded(EcalRecHitCollection::const_iterator candidate_it,
0241                                       EcalRecHitCollection::const_iterator previous_it) {
0242   // crystal should not be included...
0243   if ((candidate_it == recHits_->end()) ||                  // ...if it does not correspond to a hit
0244       (used_s.find(candidate_it->id()) != used_s.end()) ||  // ...if it already belongs to a cluster
0245       (candidate_it->energy() <= 0) ||                      // ...if it has a negative or zero energy
0246       (candidate_it->energy() > previous_it->energy()) ||   // ...or if the previous crystal had lower E
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     //      if (hit_p != 0)
0272     //  {
0273     energy += hit_p.energy();
0274     //  }
0275     //      else
0276     //  {
0277     //    std::cout << "DEBUG ALERT: Requested rechit has gone missing from rechits map! :-S" << std::endl;
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 }