From 3f213a822e55ab6a42986f18350f91b5260bafb0 Mon Sep 17 00:00:00 2001 From: Christopher Jones Date: Thu, 21 Nov 2024 08:46:26 -0600 Subject: [PATCH 1/2] Do not hold memory between events in Multi5x5ClusterAlgo - removed member variables which were just used as temporary storage - refactored variables which must be updated together to be in a class - modernized C++ syntax used --- .../interface/Multi5x5ClusterAlgo.h | 72 ++--- .../src/Multi5x5ClusterAlgo.cc | 260 ++++++++++-------- 2 files changed, 173 insertions(+), 159 deletions(-) diff --git a/RecoEcal/EgammaClusterAlgos/interface/Multi5x5ClusterAlgo.h b/RecoEcal/EgammaClusterAlgos/interface/Multi5x5ClusterAlgo.h index b1cc6dfd48b50..18af1d40265b9 100644 --- a/RecoEcal/EgammaClusterAlgos/interface/Multi5x5ClusterAlgo.h +++ b/RecoEcal/EgammaClusterAlgos/interface/Multi5x5ClusterAlgo.h @@ -21,6 +21,7 @@ #include #include #include +#include typedef std::map RecHitsMap; @@ -31,18 +32,15 @@ class Multi5x5ClusterAlgo { //so we define a proto basic cluster class which contains all the information which would be in a basic cluster //which allows the addition of its seed and the removal of a seed of another cluster easily class ProtoBasicCluster { - float energy_; - EcalRecHit seed_; std::vector > hits_; + EcalRecHit seed_; + float energy_; bool containsSeed_; public: ProtoBasicCluster(); - ProtoBasicCluster(float iEnergy, const EcalRecHit &iSeed, std::vector > &iHits) - : energy_(iEnergy), seed_(iSeed) { - hits_.swap(iHits); - containsSeed_ = isSeedCrysInHits_(); - } + ProtoBasicCluster(float iEnergy, const EcalRecHit &iSeed, std::vector > iHits) + : hits_(std::move(iHits)), seed_(iSeed), energy_(iEnergy), containsSeed_{isSeedCrysInHits_()} {} float energy() const { return energy_; } const EcalRecHit &seed() const { return seed_; } @@ -63,8 +61,8 @@ class Multi5x5ClusterAlgo { const std::vector &v_chstatus, const PositionCalc &posCalc, bool reassignSeedCrysToClusterItSeeds = false) - : ecalBarrelSeedThreshold(ebst), - ecalEndcapSeedThreshold(ecst), + : ecalBarrelSeedThreshold_(ebst), + ecalEndcapSeedThreshold_(ecst), v_chstatus_(v_chstatus), reassignSeedCrysToClusterItSeeds_(reassignSeedCrysToClusterItSeeds) { posCalculator_ = posCalc; @@ -94,60 +92,46 @@ class Multi5x5ClusterAlgo { reco::CaloID::Detectors detector_; // Energy required for a seed: - double ecalBarrelSeedThreshold; - double ecalEndcapSeedThreshold; - - // collection of all rechits - const EcalRecHitCollection *recHits_; - - // The vector of seeds: - std::vector seeds; - - std::vector > whichClusCrysBelongsTo_; - - // The set of used DetID's - std::set used_s; - std::set canSeed_s; // set of crystals not to be added but which can seed - // a new 3x3 (e.g. the outer crystals in a 5x5) - - // The vector of DetId's in the cluster currently reconstructed - std::vector > current_v; + double ecalBarrelSeedThreshold_; + double ecalEndcapSeedThreshold_; - // The vector of clusters - std::vector clusters_v; - std::vector protoClusters_; // recHit flag to be excluded from seeding std::vector v_chstatus_; bool reassignSeedCrysToClusterItSeeds_; //the seed of the 5x5 crystal is sometimes in another basic cluster, however we may want to put it back into the cluster it seeds - void mainSearch(const EcalRecHitCollection *hits, - const CaloSubdetectorGeometry *geometry_p, - const CaloSubdetectorTopology *topology_p, - const CaloSubdetectorGeometry *geometryES_p); + std::vector mainSearch(const EcalRecHitCollection *hits, + const CaloSubdetectorGeometry *geometry_p, + const CaloSubdetectorTopology *topology_p, + const CaloSubdetectorGeometry *geometryES_p, + const std::vector &seeds); // Is the crystal at the navigator position a // local maxiumum in energy? - bool checkMaxima(CaloNavigator &navigator, const EcalRecHitCollection *hits); + bool checkMaxima(CaloNavigator &navigator, const EcalRecHitCollection *hits) const; // prepare the 5x5 taking care over which crystals // are allowed to seed new clusters and which are not // after the preparation is complete - void prepareCluster(CaloNavigator &navigator, - const EcalRecHitCollection *hits, - const CaloSubdetectorGeometry *geometry); + std::vector > prepareCluster(CaloNavigator &navigator, + const EcalRecHitCollection *hits, + const CaloSubdetectorGeometry *geometry, + std::set &used_seeds, + std::set &canSeed_s) const; // Add the crystal with DetId det to the current // vector of crystals if it meets certain criteria - void addCrystal(const DetId &det); + static bool addCrystal(const DetId &det, const EcalRecHitCollection &recHits); // take the crystals in the current_v and build // them into a BasicCluster - void makeCluster(const EcalRecHitCollection *hits, - const CaloSubdetectorGeometry *geometry_p, - const CaloSubdetectorGeometry *geometryES_p, - const EcalRecHitCollection::const_iterator &seedIt, - bool seedOutside); + // NOTE: this can't be const because of posCalculator_ + std::optional makeCluster(const EcalRecHitCollection *hits, + const CaloSubdetectorGeometry *geometry_p, + const CaloSubdetectorGeometry *geometryES_p, + const EcalRecHitCollection::const_iterator &seedIt, + bool seedOutside, + std::vector > ¤t_v); }; #endif diff --git a/RecoEcal/EgammaClusterAlgos/src/Multi5x5ClusterAlgo.cc b/RecoEcal/EgammaClusterAlgos/src/Multi5x5ClusterAlgo.cc index 6480db3362c12..bce173c5978d7 100644 --- a/RecoEcal/EgammaClusterAlgos/src/Multi5x5ClusterAlgo.cc +++ b/RecoEcal/EgammaClusterAlgos/src/Multi5x5ClusterAlgo.cc @@ -58,24 +58,19 @@ std::vector Multi5x5ClusterAlgo::makeClusters(const EcalRecH reco::CaloID::Detectors detector, bool regional, const std::vector& regions) { - seeds.clear(); - used_s.clear(); - canSeed_s.clear(); - clusters_v.clear(); - - recHits_ = hits; + std::vector seeds; double threshold = 0; std::string ecalPart_string; detector_ = reco::CaloID::DET_NONE; if (detector == reco::CaloID::DET_ECAL_ENDCAP) { detector_ = reco::CaloID::DET_ECAL_ENDCAP; - threshold = ecalEndcapSeedThreshold; + threshold = ecalEndcapSeedThreshold_; ecalPart_string = "EndCap"; } if (detector == reco::CaloID::DET_ECAL_BARREL) { detector_ = reco::CaloID::DET_ECAL_BARREL; - threshold = ecalBarrelSeedThreshold; + threshold = ecalBarrelSeedThreshold_; ecalPart_string = "Barrel"; } @@ -88,20 +83,18 @@ std::vector Multi5x5ClusterAlgo::makeClusters(const EcalRecH nregions = regions.size(); if (!regional || nregions) { - EcalRecHitCollection::const_iterator it; - for (it = hits->begin(); it != hits->end(); it++) { - double energy = it->energy(); + for (auto const& hit : *hits) { + double energy = hit.energy(); if (energy < threshold) continue; // need to check to see if this line is useful! - auto thisCell = geometry_p->getGeometry(it->id()); + auto thisCell = geometry_p->getGeometry(hit.id()); // Require that RecHit is within clustering region in case // of regional reconstruction bool withinRegion = false; if (regional) { - std::vector::const_iterator region; - for (region = regions.begin(); region != regions.end(); region++) { - if (region->inRegion(thisCell->etaPos(), thisCell->phiPos())) { + for (auto const& region : regions) { + if (region.inRegion(thisCell->etaPos(), thisCell->phiPos())) { withinRegion = true; break; } @@ -109,9 +102,9 @@ std::vector Multi5x5ClusterAlgo::makeClusters(const EcalRecH } if (!regional || withinRegion) { - float ET = it->energy() * thisCell->getPosition().basicVector().unit().perp(); + float ET = hit.energy() * thisCell->getPosition().basicVector().unit().perp(); if (ET > threshold) - seeds.push_back(*it); + seeds.push_back(hit); } } } @@ -120,7 +113,7 @@ std::vector Multi5x5ClusterAlgo::makeClusters(const EcalRecH LogTrace("EcalClusters") << "Total number of seeds found in event = " << seeds.size(); - mainSearch(hits, geometry_p, topology_p, geometryES_p); + auto clusters_v = mainSearch(hits, geometry_p, topology_p, geometryES_p, seeds); sort(clusters_v.rbegin(), clusters_v.rend(), isClusterEtLess); LogTrace("EcalClusters") << "---------- end of main search. clusters have been sorted ----"; @@ -130,30 +123,91 @@ std::vector Multi5x5ClusterAlgo::makeClusters(const EcalRecH // Search for clusters // -void Multi5x5ClusterAlgo::mainSearch(const EcalRecHitCollection* hits, - const CaloSubdetectorGeometry* geometry_p, - const CaloSubdetectorTopology* topology_p, - const CaloSubdetectorGeometry* geometryES_p) { +namespace { + class ProtoClusters { + using ProtoBasicCluster = Multi5x5ClusterAlgo::ProtoBasicCluster; + std::vector clusters_; + std::vector > whichClusCrysBelongsTo_; + bool reassignSeedCrysToClusterItSeeds_; + + public: + ProtoClusters(bool reassignSeedCrysToClusterItSeeds) + : reassignSeedCrysToClusterItSeeds_(reassignSeedCrysToClusterItSeeds) {} + void push_back(ProtoBasicCluster&& c) { + clusters_.push_back(std::move(c)); + if (reassignSeedCrysToClusterItSeeds_) { + for (auto const& hit : clusters_.back().hits()) { + whichClusCrysBelongsTo_.push_back(std::pair(hit.first, clusters_.size())); + } + } + } + std::vector finalize() { + if (reassignSeedCrysToClusterItSeeds_) { + std::sort(whichClusCrysBelongsTo_.begin(), whichClusCrysBelongsTo_.end(), PairSortByFirst()); + + for (size_t clusNr = 0; clusNr < clusters_.size(); clusNr++) { + if (!clusters_[clusNr].containsSeed()) { + const EcalRecHit& seedHit = clusters_[clusNr].seed(); + typedef std::vector >::iterator It; + std::pair result = std::equal_range(whichClusCrysBelongsTo_.begin(), + whichClusCrysBelongsTo_.end(), + seedHit.id(), + PairSortByFirst()); + + if (result.first != result.second) + clusters_[result.first->second].removeHit(seedHit); + clusters_[clusNr].addSeed(); + } + } + } + whichClusCrysBelongsTo_.clear(); + return std::move(clusters_); + } + }; +} // namespace + +std::vector Multi5x5ClusterAlgo::mainSearch(const EcalRecHitCollection* hits, + const CaloSubdetectorGeometry* geometry_p, + const CaloSubdetectorTopology* topology_p, + const CaloSubdetectorGeometry* geometryES_p, + std::vector const& seeds) { LogTrace("EcalClusters") << "Building clusters............"; + // The vector of clusters + std::vector clusters_v; + + // The set of used DetID's + std::set used_s; + // set of crystals not to be added but which can seed + // a new 3x3 (e.g. the outer crystals in a 5x5) + std::set canSeed_s; + + ProtoClusters protoClusters{reassignSeedCrysToClusterItSeeds_}; + // Loop over seeds: - std::vector::iterator it; - for (it = seeds.begin(); it != seeds.end(); it++) { + bool first = true; + for (auto const& seed : seeds) { + struct Guard { + bool& b_; + Guard(bool& b) : b_{b} {}; + ~Guard() { b_ = false; } + }; + Guard guard(first); // check if this crystal is able to seed // (event though it is already used) bool usedButCanSeed = false; - if (canSeed_s.find(it->id()) != canSeed_s.end()) + if (canSeed_s.find(seed.id()) != canSeed_s.end()) usedButCanSeed = true; // avoid seeding for anomalous channels - if (!it->checkFlag(EcalRecHit::kGood)) { // if rechit is good, no need for further checks - if (it->checkFlags(v_chstatus_)) + if (!seed.checkFlag(EcalRecHit::kGood)) { // if rechit is good, no need for further checks + if (seed.checkFlags(v_chstatus_)) continue; } // make sure the current seed does not belong to a cluster already. - if ((used_s.find(it->id()) != used_s.end()) && (usedButCanSeed == false)) { - if (it == seeds.begin()) { + if ((used_s.find(seed.id()) != used_s.end()) && (usedButCanSeed == false)) { + if (first) { LogTrace("EcalClusters") << "##############################################################"; LogTrace("EcalClusters") << "DEBUG ALERT: Highest energy seed already belongs to a cluster!"; LogTrace("EcalClusters") << "##############################################################"; @@ -164,12 +218,9 @@ void Multi5x5ClusterAlgo::mainSearch(const EcalRecHitCollection* hits, continue; } - // clear the vector of hits in current cluster - current_v.clear(); - // Create a navigator at the seed and get seed // energy - CaloNavigator navigator(it->id(), topology_p); + CaloNavigator navigator(seed.id(), topology_p); DetId seedId = navigator.pos(); EcalRecHitCollection::const_iterator seedIt = hits->find(seedId); navigator.setHome(seedId); @@ -180,70 +231,57 @@ void Multi5x5ClusterAlgo::mainSearch(const EcalRecHitCollection* hits, if (localMaxima) { // build the 5x5 taking care over which crystals // can seed new clusters and which can't - prepareCluster(navigator, hits, geometry_p); - } - - // If some crystals in the current vector then - // make them into a cluster - if (!current_v.empty()) { - makeCluster(hits, geometry_p, geometryES_p, seedIt, usedButCanSeed); + auto current_v = prepareCluster(navigator, hits, geometry_p, used_s, canSeed_s); + + // If some crystals in the current vector then + // make them into a cluster + if (!current_v.empty()) { + auto c = makeCluster(hits, geometry_p, geometryES_p, seedIt, usedButCanSeed, current_v); + if (c) { + protoClusters.push_back(std::move(*c)); + } else { + for (auto const& current : current_v) { + used_s.erase(current.first); + } + } + } } - } // End loop on seed crystals - if (reassignSeedCrysToClusterItSeeds_) { - std::sort(whichClusCrysBelongsTo_.begin(), whichClusCrysBelongsTo_.end(), PairSortByFirst()); - - for (size_t clusNr = 0; clusNr < protoClusters_.size(); clusNr++) { - if (!protoClusters_[clusNr].containsSeed()) { - const EcalRecHit& seedHit = protoClusters_[clusNr].seed(); - typedef std::vector >::iterator It; - std::pair result = std::equal_range(whichClusCrysBelongsTo_.begin(), - whichClusCrysBelongsTo_.end(), - seedHit.id(), - PairSortByFirst()); - - if (result.first != result.second) - protoClusters_[result.first->second].removeHit(seedHit); - protoClusters_[clusNr].addSeed(); - } - } - } + auto finalProtoClusters = protoClusters.finalize(); - for (size_t clusNr = 0; clusNr < protoClusters_.size(); clusNr++) { - const ProtoBasicCluster& protoCluster = protoClusters_[clusNr]; - Point position; - position = posCalculator_.Calculate_Location(protoCluster.hits(), hits, geometry_p, geometryES_p); - clusters_v.push_back(reco::BasicCluster(protoCluster.energy(), - position, - reco::CaloID(detector_), - protoCluster.hits(), - reco::CaloCluster::multi5x5, - protoCluster.seed().id())); + for (auto const& protoCluster : finalProtoClusters) { + Point position = posCalculator_.Calculate_Location(protoCluster.hits(), hits, geometry_p, geometryES_p); + clusters_v.emplace_back(protoCluster.energy(), + position, + reco::CaloID(detector_), + protoCluster.hits(), + reco::CaloCluster::multi5x5, + protoCluster.seed().id()); } - protoClusters_.clear(); - whichClusCrysBelongsTo_.clear(); + return clusters_v; } -void Multi5x5ClusterAlgo::makeCluster(const EcalRecHitCollection* hits, - const CaloSubdetectorGeometry* geometry, - const CaloSubdetectorGeometry* geometryES, - const EcalRecHitCollection::const_iterator& seedIt, - bool seedOutside) { +std::optional Multi5x5ClusterAlgo::makeCluster( + const EcalRecHitCollection* hits, + const CaloSubdetectorGeometry* geometry, + const CaloSubdetectorGeometry* geometryES, + const EcalRecHitCollection::const_iterator& seedIt, + bool seedOutside, + std::vector >& current_v) { double energy = 0; //double chi2 = 0; reco::CaloID caloID; Point position; position = posCalculator_.Calculate_Location(current_v, hits, geometry, geometryES); - std::vector >::iterator it; - for (it = current_v.begin(); it != current_v.end(); it++) { - EcalRecHitCollection::const_iterator itt = hits->find((*it).first); + for (auto const& hitInfo : current_v) { + EcalRecHitCollection::const_iterator itt = hits->find(hitInfo.first); EcalRecHit hit_p = *itt; energy += hit_p.energy(); //chi2 += 0; - if ((*it).first.subdetId() == EcalBarrel) { + if (hitInfo.first.subdetId() == EcalBarrel) { caloID = reco::CaloID::DET_ECAL_BARREL; } else { caloID = reco::CaloID::DET_ECAL_ENDCAP; @@ -262,48 +300,39 @@ void Multi5x5ClusterAlgo::makeCluster(const EcalRecHitCollection* hits, // must be at least the seed energy double seedEnergy = seedIt->energy(); if ((seedOutside && energy >= 0) || (!seedOutside && energy >= seedEnergy)) { - if (reassignSeedCrysToClusterItSeeds_) { //if we're not doing this, we dont need this info so lets not bother filling it - for (size_t hitNr = 0; hitNr < current_v.size(); hitNr++) - whichClusCrysBelongsTo_.push_back(std::pair(current_v[hitNr].first, protoClusters_.size())); - } - protoClusters_.push_back(ProtoBasicCluster(energy, *seedIt, current_v)); + return ProtoBasicCluster(energy, *seedIt, std::move(current_v)); // clusters_v.push_back(reco::BasicCluster(energy, position, reco::CaloID(detector_), current_v, reco::CaloCluster::multi5x5, seedIt->id())); // if no valid cluster was built, // then free up these crystals to be used in the next... - } else { - std::vector >::iterator iter; - for (iter = current_v.begin(); iter != current_v.end(); iter++) { - used_s.erase(iter->first); - } //for(iter) } //else + return {}; } -bool Multi5x5ClusterAlgo::checkMaxima(CaloNavigator& navigator, const EcalRecHitCollection* hits) { +bool Multi5x5ClusterAlgo::checkMaxima(CaloNavigator& navigator, const EcalRecHitCollection* hits) const { bool maxima = true; EcalRecHitCollection::const_iterator thisHit; EcalRecHitCollection::const_iterator seedHit = hits->find(navigator.pos()); double seedEnergy = seedHit->energy(); - std::vector swissCrossVec; - swissCrossVec.clear(); + std::array swissCrossVec; - swissCrossVec.push_back(navigator.west()); + swissCrossVec[0] = navigator.west(); navigator.home(); - swissCrossVec.push_back(navigator.east()); + swissCrossVec[1] = navigator.east(); navigator.home(); - swissCrossVec.push_back(navigator.north()); + swissCrossVec[2] = navigator.north(); navigator.home(); - swissCrossVec.push_back(navigator.south()); + swissCrossVec[3] = navigator.south(); navigator.home(); - for (unsigned int i = 0; i < swissCrossVec.size(); ++i) { + for (auto const& det : swissCrossVec) { // look for this hit - thisHit = recHits_->find(swissCrossVec[i]); + thisHit = hits->find(det); // continue if this hit was not found - if ((swissCrossVec[i] == DetId(0)) || thisHit == recHits_->end()) + if ((det == DetId(0)) || thisHit == hits->end()) continue; // the recHit has to be skipped in the local maximum search if it was found @@ -324,9 +353,12 @@ bool Multi5x5ClusterAlgo::checkMaxima(CaloNavigator& navigator, const Eca return maxima; } -void Multi5x5ClusterAlgo::prepareCluster(CaloNavigator& navigator, - const EcalRecHitCollection* hits, - const CaloSubdetectorGeometry* geometry) { +std::vector > Multi5x5ClusterAlgo::prepareCluster(CaloNavigator& navigator, + const EcalRecHitCollection* hits, + const CaloSubdetectorGeometry* geometry, + std::set& used_s, + std::set& canSeed_s) const { + std::vector > current_v; DetId thisDet; std::set::iterator setItr; @@ -344,7 +376,10 @@ void Multi5x5ClusterAlgo::prepareCluster(CaloNavigator& navigator, // add the current crystal //std::cout << "adding " << dx << ", " << dy << std::endl; - addCrystal(thisDet); + if (addCrystal(thisDet, *hits) and (used_s.find(thisDet) == used_s.end())) { + current_v.push_back(std::pair(thisDet, 1.)); // by default hit energy fractions are set at 1. + used_s.insert(thisDet); + } // now consider if we are in an edge (outer 16) // or central (inner 9) region @@ -371,17 +406,12 @@ void Multi5x5ClusterAlgo::prepareCluster(CaloNavigator& navigator, //std::cout << "*** " << std::endl; //std::cout << " current_v contains " << current_v.size() << std::endl; //std::cout << "*** " << std::endl; + return current_v; } -void Multi5x5ClusterAlgo::addCrystal(const DetId& det) { - EcalRecHitCollection::const_iterator thisIt = recHits_->find(det); - if ((thisIt != recHits_->end()) && (thisIt->id() != DetId(0))) { - if ((used_s.find(thisIt->id()) == used_s.end())) { - //std::cout << " ... this is a good crystal and will be added" << std::endl; - current_v.push_back(std::pair(det, 1.)); // by default hit energy fractions are set at 1. - used_s.insert(det); - } - } +bool Multi5x5ClusterAlgo::addCrystal(const DetId& det, const EcalRecHitCollection& recHits) { + auto thisIt = recHits.find(det); + return ((thisIt != recHits.end()) && (thisIt->id() != DetId(0))); } bool Multi5x5ClusterAlgo::ProtoBasicCluster::removeHit(const EcalRecHit& hitToRM) { @@ -421,8 +451,8 @@ bool Multi5x5ClusterAlgo::ProtoBasicCluster::addSeed() { } bool Multi5x5ClusterAlgo::ProtoBasicCluster::isSeedCrysInHits_() const { - for (size_t hitNr = 0; hitNr < hits_.size(); hitNr++) { - if (seed_.id() == hits_[hitNr].first) + for (auto const& hit : hits_) { + if (seed_.id() == hit.first) return true; } return false; From c867babe601ef95adddfb278437f9333a46b2921 Mon Sep 17 00:00:00 2001 From: Christopher Jones Date: Thu, 21 Nov 2024 12:40:07 -0600 Subject: [PATCH 2/2] Modernize Multi5x5ClusterProducer - removed unnecessary member data - avoid unnecessary allocations - use new framework syntax - added fillDescriptions --- .../src/Multi5x5ClusterProducer.cc | 195 ++++++++---------- 1 file changed, 84 insertions(+), 111 deletions(-) diff --git a/RecoEcal/EgammaClusterProducers/src/Multi5x5ClusterProducer.cc b/RecoEcal/EgammaClusterProducers/src/Multi5x5ClusterProducer.cc index 513f227d5068d..993a5b36c6505 100644 --- a/RecoEcal/EgammaClusterProducers/src/Multi5x5ClusterProducer.cc +++ b/RecoEcal/EgammaClusterProducers/src/Multi5x5ClusterProducer.cc @@ -14,6 +14,7 @@ #include "FWCore/MessageLogger/interface/MessageLogger.h" #include "FWCore/ParameterSet/interface/ParameterSet.h" #include "FWCore/Utilities/interface/ESGetToken.h" +#include "FWCore/Utilities/interface/EDPutToken.h" #include "FWCore/Utilities/interface/Exception.h" #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h" #include "Geometry/CaloGeometry/interface/CaloGeometry.h" @@ -24,6 +25,7 @@ #include "Geometry/Records/interface/CaloGeometryRecord.h" #include "RecoEcal/EgammaClusterAlgos/interface/Multi5x5ClusterAlgo.h" #include "RecoEcal/EgammaCoreTools/interface/PositionCalc.h" +#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" #include #include @@ -34,142 +36,113 @@ class Multi5x5ClusterProducer : public edm::stream::EDProducer<> { public: Multi5x5ClusterProducer(const edm::ParameterSet& ps); - ~Multi5x5ClusterProducer() override; + void produce(edm::Event&, const edm::EventSetup&) final; - void produce(edm::Event&, const edm::EventSetup&) override; + static void fillDescriptions(edm::ConfigurationDescriptions&); private: - int nMaxPrintout_; // max # of printouts - int nEvt_; // internal counter of events + const edm::EDGetTokenT barrelHitToken_; + const edm::EDGetTokenT endcapHitToken_; + const edm::ESGetToken caloGeometryToken_; - // cluster which regions - bool doBarrel_; - bool doEndcap_; - - edm::EDGetTokenT barrelHitToken_; - edm::EDGetTokenT endcapHitToken_; - edm::ESGetToken caloGeometryToken_; - - std::string barrelClusterCollection_; - std::string endcapClusterCollection_; - - PositionCalc posCalculator_; // position calculation algorithm - Multi5x5ClusterAlgo* island_p; + const edm::EDPutTokenT barrelToken_; + const edm::EDPutTokenT endcapToken_; - bool counterExceeded() const { return ((nEvt_ > nMaxPrintout_) || (nMaxPrintout_ < 0)); } + Multi5x5ClusterAlgo island_; - const EcalRecHitCollection* getCollection(edm::Event& evt, const edm::EDGetTokenT& token); - - void clusterizeECALPart(edm::Event& evt, - const edm::EventSetup& es, - const edm::EDGetTokenT& token, - const std::string& clusterCollection, - const reco::CaloID::Detectors detector); - - void outputValidationInfo(reco::CaloClusterPtrVector& clusterPtrVector); + // cluster which regions + const bool doBarrel_; + const bool doEndcap_; + + reco::BasicClusterCollection clusterizeECALPart(const EcalRecHitCollection& hits, + const edm::EventSetup& es, + const reco::CaloID::Detectors detector); + + reco::BasicClusterCollection makeClusters(const EcalRecHitCollection& hits, + const CaloSubdetectorGeometry* geom, + const CaloSubdetectorGeometry* preshower, + const CaloSubdetectorTopology* topology, + const reco::CaloID::Detectors detector); }; #include "FWCore/Framework/interface/MakerMacros.h" DEFINE_FWK_MODULE(Multi5x5ClusterProducer); -Multi5x5ClusterProducer::Multi5x5ClusterProducer(const edm::ParameterSet& ps) { - // Parameters to identify the hit collections - barrelHitToken_ = consumes(ps.getParameter("barrelHitTag")); - - endcapHitToken_ = consumes(ps.getParameter("endcapHitTag")); - - //EventSetup Token for CaloGeometry - caloGeometryToken_ = esConsumes(); - - // should cluster algo be run in barrel and endcap? - doEndcap_ = ps.getParameter("doEndcap"); - doBarrel_ = ps.getParameter("doBarrel"); - - // The names of the produced cluster collections - barrelClusterCollection_ = ps.getParameter("barrelClusterCollection"); - endcapClusterCollection_ = ps.getParameter("endcapClusterCollection"); - - // Island algorithm parameters - double barrelSeedThreshold = ps.getParameter("IslandBarrelSeedThr"); - double endcapSeedThreshold = ps.getParameter("IslandEndcapSeedThr"); - - const std::vector flagnames = ps.getParameter >("RecHitFlagToBeExcluded"); - - const std::vector v_chstatus = StringToEnumValue(flagnames); - - // Parameters for the position calculation: - edm::ParameterSet posCalcParameters = ps.getParameter("posCalcParameters"); - posCalculator_ = PositionCalc(posCalcParameters); - - // Produces a collection of barrel and a collection of endcap clusters - produces(endcapClusterCollection_); - produces(barrelClusterCollection_); - - bool reassignSeedCrysToClusterItSeeds = false; - if (ps.exists("reassignSeedCrysToClusterItSeeds")) - reassignSeedCrysToClusterItSeeds = ps.getParameter("reassignSeedCrysToClusterItSeeds"); - - island_p = new Multi5x5ClusterAlgo( - barrelSeedThreshold, endcapSeedThreshold, v_chstatus, posCalculator_, reassignSeedCrysToClusterItSeeds); - - nEvt_ = 0; +Multi5x5ClusterProducer::Multi5x5ClusterProducer(const edm::ParameterSet& ps) + : barrelHitToken_{consumes(ps.getParameter("barrelHitTag"))}, + endcapHitToken_{consumes(ps.getParameter("endcapHitTag"))}, + caloGeometryToken_{esConsumes()}, + barrelToken_{produces(ps.getParameter("barrelClusterCollection"))}, + endcapToken_{produces(ps.getParameter("endcapClusterCollection"))}, + island_{ps.getParameter("IslandBarrelSeedThr"), + ps.getParameter("IslandEndcapSeedThr"), + StringToEnumValue(ps.getParameter>("RecHitFlagToBeExcluded")), + PositionCalc(ps.getParameter("posCalcParameters")), + ps.getParameter("reassignSeedCrysToClusterItSeeds")}, + doBarrel_{ps.getParameter("doBarrel")}, + doEndcap_{ps.getParameter("doEndcap")} {} + +void Multi5x5ClusterProducer::fillDescriptions(edm::ConfigurationDescriptions& iDesc) { + edm::ParameterSetDescription ps; + ps.add("barrelHitTag"); + ps.add("endcapHitTag"); + ps.add("doEndcap"); + ps.add("doBarrel"); + ps.add("barrelClusterCollection"); + ps.add("endcapClusterCollection"); + ps.add("IslandBarrelSeedThr"); + ps.add("IslandEndcapSeedThr"); + ps.add>("RecHitFlagToBeExcluded"); + + edm::ParameterSetDescription posCal; + posCal.add("LogWeighted"); + posCal.add("T0_barl"); + posCal.add("T0_endc"); + posCal.add("T0_endcPresh"); + posCal.add("W0"); + posCal.add("X0"); + ps.add("posCalcParameters", posCal); + + ps.add("reassignSeedCrysToClusterItSeeds", false); + iDesc.addDefault(ps); } -Multi5x5ClusterProducer::~Multi5x5ClusterProducer() { delete island_p; } - void Multi5x5ClusterProducer::produce(edm::Event& evt, const edm::EventSetup& es) { if (doEndcap_) { - clusterizeECALPart(evt, es, endcapHitToken_, endcapClusterCollection_, reco::CaloID::DET_ECAL_ENDCAP); + // get the hit collection from the event: + const EcalRecHitCollection& hitCollection = evt.get(endcapHitToken_); + evt.emplace(endcapToken_, clusterizeECALPart(hitCollection, es, reco::CaloID::DET_ECAL_ENDCAP)); } if (doBarrel_) { - clusterizeECALPart(evt, es, barrelHitToken_, barrelClusterCollection_, reco::CaloID::DET_ECAL_BARREL); + // get the hit collection from the event: + const EcalRecHitCollection& hitCollection = evt.get(barrelHitToken_); + evt.emplace(barrelToken_, clusterizeECALPart(hitCollection, es, reco::CaloID::DET_ECAL_BARREL)); } - - nEvt_++; } -const EcalRecHitCollection* Multi5x5ClusterProducer::getCollection( - edm::Event& evt, const edm::EDGetTokenT& token) { - edm::Handle rhcHandle; - evt.getByToken(token, rhcHandle); - return rhcHandle.product(); +reco::BasicClusterCollection Multi5x5ClusterProducer::makeClusters(const EcalRecHitCollection& hits, + const CaloSubdetectorGeometry* geom, + const CaloSubdetectorGeometry* preshower, + const CaloSubdetectorTopology* topology, + const reco::CaloID::Detectors detector) { + // Run the clusterization algorithm: + return island_.makeClusters(&hits, geom, topology, preshower, detector); } -void Multi5x5ClusterProducer::clusterizeECALPart(edm::Event& evt, - const edm::EventSetup& es, - const edm::EDGetTokenT& token, - const std::string& clusterCollection, - const reco::CaloID::Detectors detector) { - // get the hit collection from the event: - const EcalRecHitCollection* hitCollection_p = getCollection(evt, token); - +reco::BasicClusterCollection Multi5x5ClusterProducer::clusterizeECALPart(const EcalRecHitCollection& hitCollection, + const edm::EventSetup& es, + const reco::CaloID::Detectors detector) { // get the geometry and topology from the event setup: - edm::ESHandle geoHandle = es.getHandle(caloGeometryToken_); - - const CaloSubdetectorGeometry* geometry_p; - std::unique_ptr topology_p; + CaloGeometry const& geo = es.getData(caloGeometryToken_); + const CaloSubdetectorGeometry* preshower_p = geo.getSubdetectorGeometry(DetId::Ecal, EcalPreshower); if (detector == reco::CaloID::DET_ECAL_BARREL) { - geometry_p = geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalBarrel); - topology_p = std::make_unique(*geoHandle); + auto geometry_p = geo.getSubdetectorGeometry(DetId::Ecal, EcalBarrel); + EcalBarrelTopology topology{geo}; + return makeClusters(hitCollection, geometry_p, preshower_p, &topology, detector); } else { - geometry_p = geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalEndcap); - topology_p = std::make_unique(*geoHandle); + auto geometry_p = geo.getSubdetectorGeometry(DetId::Ecal, EcalEndcap); + EcalEndcapTopology topology{geo}; + return makeClusters(hitCollection, geometry_p, preshower_p, &topology, detector); } - - const CaloSubdetectorGeometry* geometryES_p; - geometryES_p = geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalPreshower); - - // Run the clusterization algorithm: - reco::BasicClusterCollection clusters; - clusters = island_p->makeClusters(hitCollection_p, geometry_p, topology_p.get(), geometryES_p, detector); - - // create a unique_ptr to a BasicClusterCollection, copy the barrel clusters into it and put in the Event: - auto clusters_p = std::make_unique(); - clusters_p->assign(clusters.begin(), clusters.end()); - edm::OrphanHandle bccHandle; - if (detector == reco::CaloID::DET_ECAL_BARREL) - bccHandle = evt.put(std::move(clusters_p), barrelClusterCollection_); - else - bccHandle = evt.put(std::move(clusters_p), endcapClusterCollection_); }