EgammaClusterAlgos/src/Multi5x5ClusterAlgo.cc

0001
0002 #include "RecoEcal/EgammaClusterAlgos/interface/Multi5x5ClusterAlgo.h"
0003
0004 // Geometry
0005 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0006 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0007 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0008 #include "Geometry/CaloTopology/interface/EcalEndcapTopology.h"
0009 #include "Geometry/CaloTopology/interface/EcalBarrelTopology.h"
0010 #include "RecoEcal/EgammaCoreTools/interface/ClusterEtLess.h"
0011 #include "DataFormats/CaloRecHit/interface/CaloID.h"
0012 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0013
0014 //temporary sorter, I'm sure this must exist already in CMSSW
0015 template <class T1, class T2, typename Comp = std::less<T1> >
0016 struct PairSortByFirst {
0017   Comp comp;
0018   bool operator()(const std::pair<T1, T2>& lhs, const std::pair<T1, T2>& rhs) { return comp(lhs.first, rhs.first); }
0019   bool operator()(const T1& lhs, const std::pair<T1, T2>& rhs) { return comp(lhs, rhs.first); }
0020   bool operator()(const std::pair<T1, T2>& lhs, const T1& rhs) { return comp(lhs.first, rhs); }
0021   bool operator()(const T1& lhs, const T1& rhs) { return comp(lhs, rhs); }
0022 };
0023
0024 struct EBDetIdSorter {
0025   bool operator()(const EBDetId& lhs, const EBDetId& rhs) {
0026     if (lhs.ieta() < rhs.ieta())
0027       return true;
0028     else if (lhs.ieta() > rhs.ieta())
0029       return false;
0030     else
0031       return lhs.iphi() < rhs.iphi();
0032   }
0033 };
0034
0035 struct EEDetIdSorter {
0036   bool operator()(const EEDetId& lhs, const EEDetId& rhs) {
0037     if (lhs.zside() < rhs.zside())
0038       return true;
0039     else if (lhs.zside() > rhs.zside())
0040       return false;
0041     else {  //z is equal, onto ix
0042       if (lhs.ix() < rhs.ix())
0043         return true;
0044       else if (lhs.ix() > rhs.ix())
0045         return false;
0046       else
0047         return lhs.iy() < rhs.iy();
0048     }
0049   }
0050 };
0051
0052 // Return a vector of clusters from a collection of EcalRecHits:
0053 //
0054 std::vector<reco::BasicCluster> Multi5x5ClusterAlgo::makeClusters(const EcalRecHitCollection* hits,
0055                                                                   const CaloSubdetectorGeometry* geometry_p,
0056                                                                   const CaloSubdetectorTopology* topology_p,
0057                                                                   const CaloSubdetectorGeometry* geometryES_p,
0058                                                                   reco::CaloID::Detectors detector,
0059                                                                   bool regional,
0060                                                                   const std::vector<RectangularEtaPhiRegion>& regions) {
0061   seeds.clear();
0062   used_s.clear();
0063   canSeed_s.clear();
0064   clusters_v.clear();
0065
0066   recHits_ = hits;
0067
0068   double threshold = 0;
0069   std::string ecalPart_string;
0070   detector_ = reco::CaloID::DET_NONE;
0071   if (detector == reco::CaloID::DET_ECAL_ENDCAP) {
0072     detector_ = reco::CaloID::DET_ECAL_ENDCAP;
0073     threshold = ecalEndcapSeedThreshold;
0074     ecalPart_string = "EndCap";
0075   }
0076   if (detector == reco::CaloID::DET_ECAL_BARREL) {
0077     detector_ = reco::CaloID::DET_ECAL_BARREL;
0078     threshold = ecalBarrelSeedThreshold;
0079     ecalPart_string = "Barrel";
0080   }
0081
0082   LogTrace("EcalClusters") << "-------------------------------------------------------------";
0083   LogTrace("EcalClusters") << "Island algorithm invoked for ECAL" << ecalPart_string;
0084   LogTrace("EcalClusters") << "Looking for seeds, energy threshold used = " << threshold << " GeV";
0085
0086   int nregions = 0;
0087   if (regional)
0088     nregions = regions.size();
0089
0090   if (!regional || nregions) {
0091     EcalRecHitCollection::const_iterator it;
0092     for (it = hits->begin(); it != hits->end(); it++) {
0093       double energy = it->energy();
0094       if (energy < threshold)
0095         continue;  // need to check to see if this line is useful!
0096
0097       auto thisCell = geometry_p->getGeometry(it->id());
0098       // Require that RecHit is within clustering region in case
0099       // of regional reconstruction
0100       bool withinRegion = false;
0101       if (regional) {
0102         std::vector<RectangularEtaPhiRegion>::const_iterator region;
0103         for (region = regions.begin(); region != regions.end(); region++) {
0104           if (region->inRegion(thisCell->etaPos(), thisCell->phiPos())) {
0105             withinRegion = true;
0106             break;
0107           }
0108         }
0109       }
0110
0111       if (!regional || withinRegion) {
0112         float ET = it->energy() * thisCell->getPosition().basicVector().unit().perp();
0113         if (ET > threshold)
0114           seeds.push_back(*it);
0115       }
0116     }
0117   }
0118
0119   sort(seeds.begin(), seeds.end(), [](auto const& x, auto const& y) { return (x.energy() > y.energy()); });
0120
0121   LogTrace("EcalClusters") << "Total number of seeds found in event = " << seeds.size();
0122
0123   mainSearch(hits, geometry_p, topology_p, geometryES_p);
0124   sort(clusters_v.rbegin(), clusters_v.rend(), isClusterEtLess);
0125
0126   LogTrace("EcalClusters") << "---------- end of main search. clusters have been sorted ----";
0127
0128   return clusters_v;
0129 }
0130
0131 // Search for clusters
0132 //
0133 void Multi5x5ClusterAlgo::mainSearch(const EcalRecHitCollection* hits,
0134                                      const CaloSubdetectorGeometry* geometry_p,
0135                                      const CaloSubdetectorTopology* topology_p,
0136                                      const CaloSubdetectorGeometry* geometryES_p) {
0137   LogTrace("EcalClusters") << "Building clusters............";
0138
0139   // Loop over seeds:
0140   std::vector<EcalRecHit>::iterator it;
0141   for (it = seeds.begin(); it != seeds.end(); it++) {
0142     // check if this crystal is able to seed
0143     // (event though it is already used)
0144     bool usedButCanSeed = false;
0145     if (canSeed_s.find(it->id()) != canSeed_s.end())
0146       usedButCanSeed = true;
0147
0148     // avoid seeding for anomalous channels
0149     if (!it->checkFlag(EcalRecHit::kGood)) {  // if rechit is good, no need for further checks
0150       if (it->checkFlags(v_chstatus_))
0151         continue;
0152     }
0153
0154     // make sure the current seed does not belong to a cluster already.
0155     if ((used_s.find(it->id()) != used_s.end()) && (usedButCanSeed == false)) {
0156       if (it == seeds.begin()) {
0157         LogTrace("EcalClusters") << "##############################################################";
0158         LogTrace("EcalClusters") << "DEBUG ALERT: Highest energy seed already belongs to a cluster!";
0159         LogTrace("EcalClusters") << "##############################################################";
0160       }
0161
0162       // seed crystal is used or is used and cannot seed a cluster
0163       // so continue to the next seed crystal...
0164       continue;
0165     }
0166
0167     // clear the vector of hits in current cluster
0168     current_v.clear();
0169
0170     // Create a navigator at the seed and get seed
0171     // energy
0172     CaloNavigator<DetId> navigator(it->id(), topology_p);
0173     DetId seedId = navigator.pos();
0174     EcalRecHitCollection::const_iterator seedIt = hits->find(seedId);
0175     navigator.setHome(seedId);
0176
0177     // Is the seed a local maximum?
0178     bool localMaxima = checkMaxima(navigator, hits);
0179
0180     if (localMaxima) {
0181       // build the 5x5 taking care over which crystals
0182       // can seed new clusters and which can't
0183       prepareCluster(navigator, hits, geometry_p);
0184     }
0185
0186     // If some crystals in the current vector then
0187     // make them into a cluster
0188     if (!current_v.empty()) {
0189       makeCluster(hits, geometry_p, geometryES_p, seedIt, usedButCanSeed);
0190     }
0191
0192   }  // End loop on seed crystals
0193
0194   if (reassignSeedCrysToClusterItSeeds_) {
0195     std::sort(whichClusCrysBelongsTo_.begin(), whichClusCrysBelongsTo_.end(), PairSortByFirst<DetId, int>());
0196
0197     for (size_t clusNr = 0; clusNr < protoClusters_.size(); clusNr++) {
0198       if (!protoClusters_[clusNr].containsSeed()) {
0199         const EcalRecHit& seedHit = protoClusters_[clusNr].seed();
0200         typedef std::vector<std::pair<DetId, int> >::iterator It;
0201         std::pair<It, It> result = std::equal_range(whichClusCrysBelongsTo_.begin(),
0202                                                     whichClusCrysBelongsTo_.end(),
0203                                                     seedHit.id(),
0204                                                     PairSortByFirst<DetId, int>());
0205
0206         if (result.first != result.second)
0207           protoClusters_[result.first->second].removeHit(seedHit);
0208         protoClusters_[clusNr].addSeed();
0209       }
0210     }
0211   }
0212
0213   for (size_t clusNr = 0; clusNr < protoClusters_.size(); clusNr++) {
0214     const ProtoBasicCluster& protoCluster = protoClusters_[clusNr];
0215     Point position;
0216     position = posCalculator_.Calculate_Location(protoCluster.hits(), hits, geometry_p, geometryES_p);
0217     clusters_v.push_back(reco::BasicCluster(protoCluster.energy(),
0218                                             position,
0219                                             reco::CaloID(detector_),
0220                                             protoCluster.hits(),
0221                                             reco::CaloCluster::multi5x5,
0222                                             protoCluster.seed().id()));
0223   }
0224
0225   protoClusters_.clear();
0226   whichClusCrysBelongsTo_.clear();
0227 }
0228
0229 void Multi5x5ClusterAlgo::makeCluster(const EcalRecHitCollection* hits,
0230                                       const CaloSubdetectorGeometry* geometry,
0231                                       const CaloSubdetectorGeometry* geometryES,
0232                                       const EcalRecHitCollection::const_iterator& seedIt,
0233                                       bool seedOutside) {
0234   double energy = 0;
0235   //double chi2   = 0;
0236   reco::CaloID caloID;
0237   Point position;
0238   position = posCalculator_.Calculate_Location(current_v, hits, geometry, geometryES);
0239
0240   std::vector<std::pair<DetId, float> >::iterator it;
0241   for (it = current_v.begin(); it != current_v.end(); it++) {
0242     EcalRecHitCollection::const_iterator itt = hits->find((*it).first);
0243     EcalRecHit hit_p = *itt;
0244     energy += hit_p.energy();
0245     //chi2 += 0;
0246     if ((*it).first.subdetId() == EcalBarrel) {
0247       caloID = reco::CaloID::DET_ECAL_BARREL;
0248     } else {
0249       caloID = reco::CaloID::DET_ECAL_ENDCAP;
0250     }
0251   }
0252   //chi2 /= energy;
0253
0254   LogTrace("EcalClusters") << "******** NEW CLUSTER ********";
0255   LogTrace("EcalClusters") << "No. of crystals = " << current_v.size();
0256   LogTrace("EcalClusters") << "     Energy     = " << energy;
0257   LogTrace("EcalClusters") << "     Phi        = " << position.phi();
0258   LogTrace("EcalClusters") << "     Eta " << position.eta();
0259   LogTrace("EcalClusters") << "*****************************";
0260
0261   // to be a valid cluster the cluster energy
0262   // must be at least the seed energy
0263   double seedEnergy = seedIt->energy();
0264   if ((seedOutside && energy >= 0) || (!seedOutside && energy >= seedEnergy)) {
0265     if (reassignSeedCrysToClusterItSeeds_) {  //if we're not doing this, we dont need this info so lets not bother filling it
0266       for (size_t hitNr = 0; hitNr < current_v.size(); hitNr++)
0267         whichClusCrysBelongsTo_.push_back(std::pair<DetId, int>(current_v[hitNr].first, protoClusters_.size()));
0268     }
0269     protoClusters_.push_back(ProtoBasicCluster(energy, *seedIt, current_v));
0270
0271     // clusters_v.push_back(reco::BasicCluster(energy, position, reco::CaloID(detector_), current_v, reco::CaloCluster::multi5x5, seedIt->id()));
0272
0273     // if no valid cluster was built,
0274     // then free up these crystals to be used in the next...
0275   } else {
0276     std::vector<std::pair<DetId, float> >::iterator iter;
0277     for (iter = current_v.begin(); iter != current_v.end(); iter++) {
0278       used_s.erase(iter->first);
0279     }  //for(iter)
0280   }  //else
0281 }
0282
0283 bool Multi5x5ClusterAlgo::checkMaxima(CaloNavigator<DetId>& navigator, const EcalRecHitCollection* hits) {
0284   bool maxima = true;
0285   EcalRecHitCollection::const_iterator thisHit;
0286   EcalRecHitCollection::const_iterator seedHit = hits->find(navigator.pos());
0287   double seedEnergy = seedHit->energy();
0288
0289   std::vector<DetId> swissCrossVec;
0290   swissCrossVec.clear();
0291
0292   swissCrossVec.push_back(navigator.west());
0293   navigator.home();
0294   swissCrossVec.push_back(navigator.east());
0295   navigator.home();
0296   swissCrossVec.push_back(navigator.north());
0297   navigator.home();
0298   swissCrossVec.push_back(navigator.south());
0299   navigator.home();
0300
0301   for (unsigned int i = 0; i < swissCrossVec.size(); ++i) {
0302     // look for this hit
0303     thisHit = recHits_->find(swissCrossVec[i]);
0304
0305     // continue if this hit was not found
0306     if ((swissCrossVec[i] == DetId(0)) || thisHit == recHits_->end())
0307       continue;
0308
0309     // the recHit has to be skipped in the local maximum search if it was found
0310     // in the map of channels to be excluded
0311     uint32_t rhFlag = thisHit->recoFlag();
0312     std::vector<int>::const_iterator vit = std::find(v_chstatus_.begin(), v_chstatus_.end(), rhFlag);
0313     if (vit != v_chstatus_.end())
0314       continue;
0315
0316     // if this crystal has more energy than the seed then we do
0317     // not have a local maxima
0318     if (thisHit->energy() > seedEnergy) {
0319       maxima = false;
0320       break;
0321     }
0322   }
0323
0324   return maxima;
0325 }
0326
0327 void Multi5x5ClusterAlgo::prepareCluster(CaloNavigator<DetId>& navigator,
0328                                          const EcalRecHitCollection* hits,
0329                                          const CaloSubdetectorGeometry* geometry) {
0330   DetId thisDet;
0331   std::set<DetId>::iterator setItr;
0332
0333   // now add the 5x5 taking care to mark the edges
0334   // as able to seed and where overlapping in the central
0335   // region with crystals that were previously able to seed
0336   // change their status so they are not able to seed
0337   //std::cout << std::endl;
0338   for (int dx = -2; dx < 3; ++dx) {
0339     for (int dy = -2; dy < 3; ++dy) {
0340       // navigate in free steps forming
0341       // a full 5x5
0342       thisDet = navigator.offsetBy(dx, dy);
0343       navigator.home();
0344
0345       // add the current crystal
0346       //std::cout << "adding " << dx << ", " << dy << std::endl;
0347       addCrystal(thisDet);
0348
0349       // now consider if we are in an edge (outer 16)
0350       // or central (inner 9) region
0351       if ((abs(dx) > 1) || (abs(dy) > 1)) {
0352         // this is an "edge" so should be allowed to seed
0353         // provided it is not already used
0354         //std::cout << "   setting can seed" << std::endl;
0355         canSeed_s.insert(thisDet);
0356       }  // end if "edge"
0357       else {
0358         // or else we are in the central 3x3
0359         // and must remove any of these crystals from the canSeed set
0360         setItr = canSeed_s.find(thisDet);
0361         if (setItr != canSeed_s.end()) {
0362           //std::cout << "   unsetting can seed" << std::endl;
0363           canSeed_s.erase(setItr);
0364         }
0365       }  // end if "centre"
0366
0367     }  // end loop on dy
0368
0369   }  // end loop on dx
0370
0371   //std::cout << "*** " << std::endl;
0372   //std::cout << " current_v contains " << current_v.size() << std::endl;
0373   //std::cout << "*** " << std::endl;
0374 }
0375
0376 void Multi5x5ClusterAlgo::addCrystal(const DetId& det) {
0377   EcalRecHitCollection::const_iterator thisIt = recHits_->find(det);
0378   if ((thisIt != recHits_->end()) && (thisIt->id() != DetId(0))) {
0379     if ((used_s.find(thisIt->id()) == used_s.end())) {
0380       //std::cout << "   ... this is a good crystal and will be added" << std::endl;
0381       current_v.push_back(std::pair<DetId, float>(det, 1.));  // by default hit energy fractions are set at 1.
0382       used_s.insert(det);
0383     }
0384   }
0385 }
0386
0387 bool Multi5x5ClusterAlgo::ProtoBasicCluster::removeHit(const EcalRecHit& hitToRM) {
0388   std::vector<std::pair<DetId, float> >::iterator hitPos;
0389   for (hitPos = hits_.begin(); hitPos < hits_.end(); hitPos++) {
0390     if (hitToRM.id() == hitPos->first)
0391       break;
0392   }
0393   if (hitPos != hits_.end()) {
0394     hits_.erase(hitPos);
0395     energy_ -= hitToRM.energy();
0396     return true;
0397   }
0398   return false;
0399 }
0400
0401 //now the tricky thing is to make sure we insert it in the right place in the hit vector
0402 //order is from -2,-2, -2,-1 to 2,2 with seed being 0,0 (eta,phi)
0403 bool Multi5x5ClusterAlgo::ProtoBasicCluster::addSeed() {
0404   typedef std::vector<std::pair<DetId, float> >::iterator It;
0405   std::pair<It, It> hitPos;
0406
0407   if (seed_.id().subdetId() == EcalBarrel) {
0408     hitPos = std::equal_range(hits_.begin(), hits_.end(), seed_.id(), PairSortByFirst<DetId, float, EBDetIdSorter>());
0409   } else {
0410     hitPos = std::equal_range(hits_.begin(), hits_.end(), seed_.id(), PairSortByFirst<DetId, float, EEDetIdSorter>());
0411   }
0412
0413   if (hitPos.first == hitPos.second) {  //it doesnt already exist in the vec, add it
0414     hits_.insert(hitPos.first, std::pair<DetId, float>(seed_.id(), 1.));
0415     energy_ += seed_.energy();
0416     containsSeed_ = true;
0417
0418     return true;
0419   } else
0420     return false;
0421 }
0422
0423 bool Multi5x5ClusterAlgo::ProtoBasicCluster::isSeedCrysInHits_() const {
0424   for (size_t hitNr = 0; hitNr < hits_.size(); hitNr++) {
0425     if (seed_.id() == hits_[hitNr].first)
0426       return true;
0427   }
0428   return false;
0429 }