Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:24:39 

0001 #include "RecoEcal/EgammaClusterAlgos/interface/PFECALSuperClusterAlgo.h"
0002 #include "RecoEcal/EgammaCoreTools/interface/EcalClustersGraph.h"
0003 #include "CommonTools/ParticleFlow/interface/PFClusterWidthAlgo.h"
0004 #include "RecoParticleFlow/PFClusterTools/interface/LinkByRecHit.h"
0005 #include "DataFormats/ParticleFlowReco/interface/PFLayer.h"
0006 #include "DataFormats/ParticleFlowReco/interface/PFRecHit.h"
0007 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0008 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0009 #include "DataFormats/EcalDetId/interface/ESDetId.h"
0010 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0011 #include "RecoEcal/EgammaCoreTools/interface/Mustache.h"
0012 #include "Math/GenVector/VectorUtil.h"
0013 #include "TVector2.h"
0014 
0015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0016 
0017 #include <memory>
0018 
0019 #include <cmath>
0020 #include <functional>
0021 #include <sstream>
0022 #include <stdexcept>
0023 #include <string>
0024 
0025 using namespace std;
0026 using namespace std::placeholders;  // for _1, _2, _3...
0027 
0028 namespace {
0029   typedef edm::View<reco::PFCluster> PFClusterView;
0030   typedef edm::Ptr<reco::PFCluster> PFClusterPtr;
0031   typedef edm::PtrVector<reco::PFCluster> PFClusterPtrVector;
0032   typedef std::pair<reco::CaloClusterPtr::key_type, reco::CaloClusterPtr> EEPSPair;
0033 
0034   bool sortByKey(const EEPSPair& a, const EEPSPair& b) { return a.first < b.first; }
0035 
0036   inline double ptFast(const double energy, const math::XYZPoint& position, const math::XYZPoint& origin) {
0037     const auto v = position - origin;
0038     return energy * std::sqrt(v.perp2() / v.mag2());
0039   }
0040 
0041   bool greaterByEt(const CalibratedPFCluster& x, const CalibratedPFCluster& y) {
0042     const math::XYZPoint zero(0, 0, 0);
0043     const double xpt = ptFast(x.energy(), x.ptr()->position(), zero);
0044     const double ypt = ptFast(y.energy(), y.ptr()->position(), zero);
0045     return xpt > ypt;
0046   }
0047 
0048   bool isSeed(const CalibratedPFCluster& x, double threshold, bool useETcut) {
0049     const math::XYZPoint zero(0, 0, 0);
0050     double e_or_et = x.energy();
0051     if (useETcut)
0052       e_or_et = ptFast(e_or_et, x.ptr()->position(), zero);
0053     return e_or_et > threshold;
0054   }
0055 
0056   bool isLinkedByRecHit(const CalibratedPFCluster& x,
0057                         const CalibratedPFCluster& seed,
0058                         const double threshold,
0059                         const double majority,
0060                         const double maxDEta,
0061                         const double maxDPhi) {
0062     if (seed.energy_nocalib() < threshold) {
0063       return false;
0064     }
0065     const double dEta = std::abs(seed.eta() - x.eta());
0066     const double dPhi = std::abs(TVector2::Phi_mpi_pi(seed.phi() - x.phi()));
0067     if (maxDEta < dEta || maxDPhi < dPhi) {
0068       return false;
0069     }
0070     // now see if the clusters overlap in rechits
0071     const auto& seedHitsAndFractions = seed.ptr()->hitsAndFractions();
0072     const auto& xHitsAndFractions = x.ptr()->hitsAndFractions();
0073     double x_rechits_tot = xHitsAndFractions.size();
0074     double x_rechits_match = 0.0;
0075     for (const std::pair<DetId, float>& seedHit : seedHitsAndFractions) {
0076       for (const std::pair<DetId, float>& xHit : xHitsAndFractions) {
0077         if (seedHit.first == xHit.first) {
0078           x_rechits_match += 1.0;
0079         }
0080       }
0081     }
0082     return x_rechits_match / x_rechits_tot > majority;
0083   }
0084 
0085   bool isClustered(const CalibratedPFCluster& x,
0086                    const CalibratedPFCluster seed,
0087                    const PFECALSuperClusterAlgo::clustering_type type,
0088                    const EcalMustacheSCParameters* mustache_params,
0089                    const EcalSCDynamicDPhiParameters* dynamic_dphi_params,
0090                    const bool dyn_dphi,
0091                    const double etawidthSuperCluster,
0092                    const double phiwidthSuperCluster) {
0093     const double dphi = std::abs(TVector2::Phi_mpi_pi(seed.phi() - x.phi()));
0094     const bool passes_dphi =
0095         ((!dyn_dphi && dphi < phiwidthSuperCluster) ||
0096          (dyn_dphi && reco::MustacheKernel::inDynamicDPhiWindow(
0097                           dynamic_dphi_params, seed.eta(), seed.phi(), x.energy_nocalib(), x.eta(), x.phi())));
0098 
0099     if (type == PFECALSuperClusterAlgo::kBOX) {
0100       return (std::abs(seed.eta() - x.eta()) < etawidthSuperCluster && passes_dphi);
0101     }
0102     if (type == PFECALSuperClusterAlgo::kMustache) {
0103       return (passes_dphi && reco::MustacheKernel::inMustache(
0104                                  mustache_params, seed.eta(), seed.phi(), x.energy_nocalib(), x.eta(), x.phi()));
0105     }
0106     return false;
0107   }
0108 
0109 }  // namespace
0110 
0111 PFECALSuperClusterAlgo::PFECALSuperClusterAlgo(const reco::SCProducerCache* cache)
0112     : beamSpot_(nullptr), SCProducerCache_(cache) {}
0113 
0114 void PFECALSuperClusterAlgo::setPFClusterCalibration(const std::shared_ptr<PFEnergyCalibration>& calib) {
0115   _pfEnergyCalibration = calib;
0116 }
0117 
0118 void PFECALSuperClusterAlgo::setTokens(const edm::ParameterSet& iConfig, edm::ConsumesCollector&& cc) {
0119   inputTagPFClusters_ = cc.consumes<edm::View<reco::PFCluster>>(iConfig.getParameter<edm::InputTag>("PFClusters"));
0120   inputTagPFClustersES_ =
0121       cc.consumes<reco::PFCluster::EEtoPSAssociation>(iConfig.getParameter<edm::InputTag>("ESAssociation"));
0122   inputTagBeamSpot_ = cc.consumes<reco::BeamSpot>(iConfig.getParameter<edm::InputTag>("BeamSpot"));
0123 
0124   esEEInterCalibToken_ =
0125       cc.esConsumes<ESEEIntercalibConstants, ESEEIntercalibConstantsRcd, edm::Transition::BeginLuminosityBlock>();
0126   esChannelStatusToken_ = cc.esConsumes<ESChannelStatus, ESChannelStatusRcd, edm::Transition::BeginLuminosityBlock>();
0127 
0128   if (_clustype == PFECALSuperClusterAlgo::kMustache) {
0129     ecalMustacheSCParametersToken_ = cc.esConsumes<EcalMustacheSCParameters, EcalMustacheSCParametersRcd>();
0130   }
0131   if (useDynamicDPhi_) {
0132     ecalSCDynamicDPhiParametersToken_ = cc.esConsumes<EcalSCDynamicDPhiParameters, EcalSCDynamicDPhiParametersRcd>();
0133   }
0134 
0135   if (useRegression_) {
0136     const edm::ParameterSet& regconf = iConfig.getParameter<edm::ParameterSet>("regressionConfig");
0137 
0138     regr_ = std::make_unique<SCEnergyCorrectorSemiParm>();
0139     regr_->setTokens(regconf, cc);
0140   }
0141 
0142   if (isOOTCollection_ || _clustype == PFECALSuperClusterAlgo::kDeepSC) {  // OOT photons and  DeepSC uses rechits
0143     inputTagBarrelRecHits_ = cc.consumes<EcalRecHitCollection>(iConfig.getParameter<edm::InputTag>("barrelRecHits"));
0144     inputTagEndcapRecHits_ = cc.consumes<EcalRecHitCollection>(iConfig.getParameter<edm::InputTag>("endcapRecHits"));
0145   }
0146   if (_clustype == PFECALSuperClusterAlgo::kDeepSC) {  // DeepSC uses geometry also
0147     caloTopologyToken_ = cc.esConsumes<CaloTopology, CaloTopologyRecord, edm::Transition::BeginLuminosityBlock>();
0148     caloGeometryToken_ = cc.esConsumes<CaloGeometry, CaloGeometryRecord, edm::Transition::BeginLuminosityBlock>();
0149   }
0150 }
0151 
0152 void PFECALSuperClusterAlgo::update(const edm::EventSetup& setup) {
0153   if (useRegression_) {
0154     regr_->setEventSetup(setup);
0155   }
0156 
0157   edm::ESHandle<ESEEIntercalibConstants> esEEInterCalibHandle_ = setup.getHandle(esEEInterCalibToken_);
0158   _pfEnergyCalibration->initAlphaGamma_ESplanes_fromDB(esEEInterCalibHandle_.product());
0159 
0160   edm::ESHandle<ESChannelStatus> esChannelStatusHandle_ = setup.getHandle(esChannelStatusToken_);
0161   channelStatus_ = esChannelStatusHandle_.product();
0162 
0163   if (_clustype == PFECALSuperClusterAlgo::kDeepSC) {  // DeepSC uses geometry
0164     edm::ESHandle<CaloGeometry> caloGeometryHandle_ = setup.getHandle(caloGeometryToken_);
0165     geometry_ = caloGeometryHandle_.product();
0166     ebGeom_ = caloGeometryHandle_->getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
0167     eeGeom_ = caloGeometryHandle_->getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
0168     esGeom_ = caloGeometryHandle_->getSubdetectorGeometry(DetId::Ecal, EcalPreshower);
0169 
0170     edm::ESHandle<CaloTopology> caloTopologyHandle_ = setup.getHandle(caloTopologyToken_);
0171     topology_ = caloTopologyHandle_.product();
0172   }
0173 }
0174 
0175 void PFECALSuperClusterAlgo::updateSCParams(const edm::EventSetup& setup) {
0176   if (_clustype == PFECALSuperClusterAlgo::kMustache) {
0177     mustacheSCParams_ = &setup.getData(ecalMustacheSCParametersToken_);
0178   }
0179   if (useDynamicDPhi_) {
0180     scDynamicDPhiParams_ = &setup.getData(ecalSCDynamicDPhiParametersToken_);
0181   }
0182 }
0183 
0184 void PFECALSuperClusterAlgo::loadAndSortPFClusters(const edm::Event& iEvent) {
0185   //load input collections
0186   //Load the pfcluster collections
0187   edm::Handle<edm::View<reco::PFCluster>> pfclustersHandle;
0188   iEvent.getByToken(inputTagPFClusters_, pfclustersHandle);
0189 
0190   edm::Handle<reco::PFCluster::EEtoPSAssociation> psAssociationHandle;
0191   iEvent.getByToken(inputTagPFClustersES_, psAssociationHandle);
0192 
0193   const PFClusterView& clusters = *pfclustersHandle.product();
0194   const reco::PFCluster::EEtoPSAssociation& psclusters = *psAssociationHandle.product();
0195 
0196   //load BeamSpot
0197   edm::Handle<reco::BeamSpot> bsHandle;
0198   iEvent.getByToken(inputTagBeamSpot_, bsHandle);
0199   beamSpot_ = bsHandle.product();
0200 
0201   //initialize regression for this event
0202   if (useRegression_) {
0203     regr_->setEvent(iEvent);
0204   }
0205 
0206   // reset the system for running
0207   superClustersEB_ = std::make_unique<reco::SuperClusterCollection>();
0208   _clustersEB.clear();
0209   superClustersEE_ = std::make_unique<reco::SuperClusterCollection>();
0210   _clustersEE.clear();
0211   EEtoPS_ = &psclusters;
0212 
0213   //Select PF clusters available for the clustering
0214   for (size_t i = 0; i < clusters.size(); ++i) {
0215     auto cluster = clusters.ptrAt(i);
0216     LogDebug("PFClustering") << "Loading PFCluster i=" << cluster.key() << " energy=" << cluster->energy();
0217 
0218     // protection for sim clusters
0219     if (cluster->caloID().detectors() == 0 && cluster->hitsAndFractions().empty())
0220       continue;
0221 
0222     CalibratedPFCluster calib_cluster(cluster);
0223     switch (cluster->layer()) {
0224       case PFLayer::ECAL_BARREL:
0225         if (calib_cluster.energy() > threshPFClusterBarrel_) {
0226           _clustersEB.push_back(calib_cluster);
0227         }
0228         break;
0229       case PFLayer::HGCAL:
0230       case PFLayer::ECAL_ENDCAP:
0231         if (calib_cluster.energy() > threshPFClusterEndcap_) {
0232           _clustersEE.push_back(calib_cluster);
0233         }
0234         break;
0235       default:
0236         break;
0237     }
0238   }
0239   // sort full cluster collections by their calibrated energy
0240   // this will put all the seeds first by construction
0241   std::sort(_clustersEB.begin(), _clustersEB.end(), greaterByEt);
0242   std::sort(_clustersEE.begin(), _clustersEE.end(), greaterByEt);
0243 
0244   // set recHit collections for OOT photons and DeepSC
0245   if (isOOTCollection_ || _clustype == PFECALSuperClusterAlgo::kDeepSC) {
0246     edm::Handle<EcalRecHitCollection> barrelRecHitsHandle;
0247     iEvent.getByToken(inputTagBarrelRecHits_, barrelRecHitsHandle);
0248     if (!barrelRecHitsHandle.isValid()) {
0249       throw cms::Exception("PFECALSuperClusterAlgo")
0250           << "If you use OOT photons or DeepSC, need to specify proper barrel rec hit collection";
0251     }
0252     barrelRecHits_ = barrelRecHitsHandle.product();
0253 
0254     edm::Handle<EcalRecHitCollection> endcapRecHitsHandle;
0255     iEvent.getByToken(inputTagEndcapRecHits_, endcapRecHitsHandle);
0256     if (!endcapRecHitsHandle.isValid()) {
0257       throw cms::Exception("PFECALSuperClusterAlgo")
0258           << "If you use OOT photons or DeepSC, need to specify proper endcap rec hit collection";
0259     }
0260     endcapRecHits_ = endcapRecHitsHandle.product();
0261   }
0262 }
0263 
0264 void PFECALSuperClusterAlgo::run() {
0265   // clusterize the EB
0266   buildAllSuperClusters(_clustersEB, threshPFClusterSeedBarrel_);
0267   // clusterize the EE
0268   buildAllSuperClusters(_clustersEE, threshPFClusterSeedEndcap_);
0269 }
0270 
0271 void PFECALSuperClusterAlgo::buildAllSuperClusters(CalibratedPFClusterVector& clusters, double seedthresh) {
0272   if (_clustype == PFECALSuperClusterAlgo::kMustache || _clustype == PFECALSuperClusterAlgo::kBOX)
0273     buildAllSuperClustersMustacheOrBox(clusters, seedthresh);
0274   else if (_clustype == PFECALSuperClusterAlgo::kDeepSC)
0275     buildAllSuperClustersDeepSC(clusters, seedthresh);
0276 }
0277 
0278 void PFECALSuperClusterAlgo::buildAllSuperClustersMustacheOrBox(CalibratedPFClusterVector& clusters,
0279                                                                 double seedthresh) {
0280   auto seedable = std::bind(isSeed, _1, seedthresh, threshIsET_);
0281 
0282   // make sure only seeds appear at the front of the list of clusters
0283   std::stable_partition(clusters.begin(), clusters.end(), seedable);
0284 
0285   // in each iteration we are working on a list that is already sorted
0286   // in the cluster energy and remains so through each iteration
0287   // NB: since clusters is sorted in loadClusters any_of has O(1)
0288   //     timing until you run out of seeds!
0289   while (std::any_of(clusters.cbegin(), clusters.cend(), seedable)) {
0290     buildSuperClusterMustacheOrBox(clusters.front(), clusters);
0291   }
0292 }
0293 
0294 void PFECALSuperClusterAlgo::buildAllSuperClustersDeepSC(CalibratedPFClusterVector& clusters, double seedthresh) {
0295   auto seedable = std::bind(isSeed, _1, seedthresh, threshIsET_);
0296   // EcalClustersGraph utility class for DeepSC algorithm application
0297   // make sure only seeds appear at the front of the list of clusters
0298   auto last_seed = std::stable_partition(clusters.begin(), clusters.end(), seedable);
0299 
0300   reco::EcalClustersGraph ecalClusterGraph_{clusters,
0301                                             static_cast<int>(std::distance(clusters.begin(), last_seed)),
0302                                             topology_,
0303                                             ebGeom_,
0304                                             eeGeom_,
0305                                             barrelRecHits_,
0306                                             endcapRecHits_,
0307                                             SCProducerCache_};
0308   // Build sub-regions of the detector where the DeepSC algo will be run
0309   ecalClusterGraph_.initWindows();
0310   // For each sub-region, prepare the DeepSC input tensors
0311   ecalClusterGraph_.fillVariables();
0312   // Evaluate the DeepSC algorithm and save the scores
0313   ecalClusterGraph_.evaluateScores();
0314   // Select the final SuperCluster using the CollectionStrategy defined in the cfi
0315   ecalClusterGraph_.setThresholds();
0316   ecalClusterGraph_.selectClusters();
0317   // Extract the final SuperCluster collection
0318   reco::EcalClustersGraph::EcalGraphOutput windows = ecalClusterGraph_.getGraphOutput();
0319   for (auto& [seed, clustered] : windows) {
0320     bool isEE = false;
0321     if (seed.ptr()->layer() == PFLayer::ECAL_ENDCAP)
0322       isEE = true;
0323     finalizeSuperCluster(seed, clustered, isEE);
0324   }
0325 }
0326 
0327 void PFECALSuperClusterAlgo::buildSuperClusterMustacheOrBox(CalibratedPFCluster& seed,
0328                                                             CalibratedPFClusterVector& clusters) {
0329   CalibratedPFClusterVector clustered;
0330 
0331   double etawidthSuperCluster = 0.0;
0332   double phiwidthSuperCluster = 0.0;
0333   bool isEE = false;
0334   switch (seed.ptr()->layer()) {
0335     case PFLayer::ECAL_BARREL:
0336       phiwidthSuperCluster = phiwidthSuperClusterBarrel_;
0337       etawidthSuperCluster = etawidthSuperClusterBarrel_;
0338       edm::LogInfo("PFClustering") << "Building SC number " << superClustersEB_->size() + 1 << " in the ECAL barrel!";
0339       break;
0340     case PFLayer::HGCAL:
0341     case PFLayer::ECAL_ENDCAP:
0342 
0343       phiwidthSuperCluster = phiwidthSuperClusterEndcap_;
0344       etawidthSuperCluster = etawidthSuperClusterEndcap_;
0345       edm::LogInfo("PFClustering") << "Building SC number " << superClustersEE_->size() + 1 << " in the ECAL endcap!"
0346                                    << std::endl;
0347       isEE = true;
0348       break;
0349     default:
0350       break;
0351   }
0352 
0353   auto isClusteredWithSeed = std::bind(isClustered,
0354                                        _1,
0355                                        seed,
0356                                        _clustype,
0357                                        mustacheSCParams_,
0358                                        scDynamicDPhiParams_,
0359                                        useDynamicDPhi_,
0360                                        etawidthSuperCluster,
0361                                        phiwidthSuperCluster);
0362 
0363   auto matchesSeedByRecHit = std::bind(isLinkedByRecHit, _1, seed, satelliteThreshold_, fractionForMajority_, 0.1, 0.2);
0364 
0365   // this function shuffles the list of clusters into a list
0366   // where all clustered sub-clusters are at the front
0367   // and returns a pointer to the first unclustered cluster.
0368   // The relative ordering of clusters is preserved
0369   // (i.e. both resulting sub-lists are sorted by energy).
0370   auto not_clustered = std::stable_partition(clusters.begin(), clusters.end(), isClusteredWithSeed);
0371   // satellite cluster merging
0372   // it was found that large clusters can split!
0373   if (doSatelliteClusterMerge_) {
0374     not_clustered = std::stable_partition(not_clustered, clusters.end(), matchesSeedByRecHit);
0375   }
0376 
0377   if (verbose_) {
0378     edm::LogInfo("PFClustering") << "Dumping cluster detail";
0379     edm::LogVerbatim("PFClustering") << "\tPassed seed: e = " << seed.energy_nocalib() << " eta = " << seed.eta()
0380                                      << " phi = " << seed.phi() << std::endl;
0381     for (auto clus = clusters.cbegin(); clus != not_clustered; ++clus) {
0382       edm::LogVerbatim("PFClustering") << "\t\tClustered cluster: e = " << (*clus).energy_nocalib()
0383                                        << " eta = " << (*clus).eta() << " phi = " << (*clus).phi() << std::endl;
0384     }
0385     for (auto clus = not_clustered; clus != clusters.end(); ++clus) {
0386       edm::LogVerbatim("PFClustering") << "\tNon-Clustered cluster: e = " << (*clus).energy_nocalib()
0387                                        << " eta = " << (*clus).eta() << " phi = " << (*clus).phi() << std::endl;
0388     }
0389   }
0390 
0391   if (not_clustered == clusters.begin()) {
0392     if (dropUnseedable_) {
0393       clusters.erase(clusters.begin());
0394       return;
0395     } else {
0396       throw cms::Exception("PFECALSuperClusterAlgo::buildSuperCluster")
0397           << "Cluster is not seedable!" << std::endl
0398           << "\tNon-Clustered cluster: e = " << (*not_clustered).energy_nocalib() << " eta = " << (*not_clustered).eta()
0399           << " phi = " << (*not_clustered).phi() << std::endl;
0400     }
0401   }
0402 
0403   // move the clustered clusters out of available cluster list
0404   // and into a temporary vector for building the SC
0405   CalibratedPFClusterVector clustered_tmp(clusters.begin(), not_clustered);
0406   clustered = clustered_tmp;
0407   clusters.erase(clusters.begin(), not_clustered);
0408 
0409   // Finalize the SuperCluster passing the list of clustered clusters
0410   finalizeSuperCluster(seed, clustered, isEE);
0411 }
0412 
0413 void PFECALSuperClusterAlgo::finalizeSuperCluster(CalibratedPFCluster& seed,
0414                                                   CalibratedPFClusterVector& clustered,
0415                                                   bool isEE) {
0416   // need the vector of raw pointers for a PF width class
0417   std::vector<const reco::PFCluster*> bare_ptrs;
0418   // calculate necessary parameters and build the SC
0419   double posX(0), posY(0), posZ(0), corrSCEnergy(0), corrPS1Energy(0), corrPS2Energy(0), energyweight(0),
0420       energyweighttot(0);
0421   for (const auto& clus : clustered) {
0422     double ePS1 = 0.0;
0423     double ePS2 = 0.0;
0424     energyweight = clus.energy_nocalib();
0425     bare_ptrs.push_back(clus.ptr().get());
0426     // update EE calibrated super cluster energies
0427     if (isEE) {
0428       auto ee_key_val = std::make_pair(clus.ptr().key(), edm::Ptr<reco::PFCluster>());
0429       const auto clustops = std::equal_range(EEtoPS_->begin(), EEtoPS_->end(), ee_key_val, sortByKey);
0430       std::vector<reco::PFCluster const*> psClusterPointers;
0431       for (auto i_ps = clustops.first; i_ps != clustops.second; ++i_ps) {
0432         psClusterPointers.push_back(i_ps->second.get());
0433       }
0434       auto calibratedEnergies = _pfEnergyCalibration->calibrateEndcapClusterEnergies(
0435           *(clus.ptr()), psClusterPointers, *channelStatus_, applyCrackCorrections_);
0436       ePS1 = calibratedEnergies.ps1Energy;
0437       ePS2 = calibratedEnergies.ps2Energy;
0438     }
0439 
0440     if (ePS1 == -1.)
0441       ePS1 = 0;
0442     if (ePS2 == -1.)
0443       ePS2 = 0;
0444 
0445     switch (_eweight) {
0446       case kRaw:  // energyweight is initialized to raw cluster energy
0447         break;
0448       case kCalibratedNoPS:
0449         energyweight = clus.energy() - ePS1 - ePS2;
0450         break;
0451       case kCalibratedTotal:
0452         energyweight = clus.energy();
0453         break;
0454       default:
0455         break;
0456     }
0457     const math::XYZPoint& cluspos = clus.ptr()->position();
0458     posX += energyweight * cluspos.X();
0459     posY += energyweight * cluspos.Y();
0460     posZ += energyweight * cluspos.Z();
0461 
0462     energyweighttot += energyweight;
0463     corrSCEnergy += clus.energy();
0464     corrPS1Energy += ePS1;
0465     corrPS2Energy += ePS2;
0466   }
0467   posX /= energyweighttot;
0468   posY /= energyweighttot;
0469   posZ /= energyweighttot;
0470 
0471   // now build the supercluster
0472   reco::SuperCluster new_sc(corrSCEnergy, math::XYZPoint(posX, posY, posZ));
0473   new_sc.setCorrectedEnergy(corrSCEnergy);
0474   new_sc.setSeed(clustered.front().ptr());
0475   new_sc.setPreshowerEnergy(corrPS1Energy + corrPS2Energy);
0476   new_sc.setPreshowerEnergyPlane1(corrPS1Energy);
0477   new_sc.setPreshowerEnergyPlane2(corrPS2Energy);
0478   for (const auto& clus : clustered) {
0479     new_sc.addCluster(clus.ptr());
0480 
0481     auto& hits_and_fractions = clus.ptr()->hitsAndFractions();
0482     for (auto& hit_and_fraction : hits_and_fractions) {
0483       new_sc.addHitAndFraction(hit_and_fraction.first, hit_and_fraction.second);
0484     }
0485     if (isEE) {
0486       auto ee_key_val = std::make_pair(clus.ptr().key(), edm::Ptr<reco::PFCluster>());
0487       const auto clustops = std::equal_range(EEtoPS_->begin(), EEtoPS_->end(), ee_key_val, sortByKey);
0488       // EE rechits should be uniquely matched to sets of pre-shower
0489       // clusters at this point, so we throw an exception if otherwise
0490       // now wrapped in EDM debug flags
0491       for (auto i_ps = clustops.first; i_ps != clustops.second; ++i_ps) {
0492         edm::Ptr<reco::PFCluster> psclus(i_ps->second);
0493 #ifdef EDM_ML_DEBUG
0494 
0495         auto found_pscluster =
0496             std::find(new_sc.preshowerClustersBegin(), new_sc.preshowerClustersEnd(), reco::CaloClusterPtr(psclus));
0497 
0498         if (found_pscluster == new_sc.preshowerClustersEnd()) {
0499 #endif
0500           new_sc.addPreshowerCluster(psclus);
0501 #ifdef EDM_ML_DEBUG
0502         } else {
0503           throw cms::Exception("PFECALSuperClusterAlgo::buildSuperCluster")
0504               << "Found a PS cluster matched to more than one EE cluster!" << std::endl
0505               << std::hex << psclus.get() << " == " << found_pscluster->get() << std::dec << std::endl;
0506         }
0507 #endif
0508       }
0509     }
0510   }
0511 
0512   // calculate linearly weighted cluster widths
0513   PFClusterWidthAlgo pfwidth(bare_ptrs);
0514   new_sc.setEtaWidth(pfwidth.pflowEtaWidth());
0515   new_sc.setPhiWidth(pfwidth.pflowPhiWidth());
0516 
0517   // cache the value of the raw energy
0518   new_sc.rawEnergy();
0519 
0520   //apply regression energy corrections
0521   if (useRegression_) {
0522     regr_->modifyObject(new_sc);
0523   }
0524 
0525   // save the super cluster to the appropriate list (if it passes the final
0526   // Et threshold)
0527   //Note that Et is computed here with respect to the beamspot position
0528   //in order to be consistent with the cut applied in the
0529   //ElectronSeedProducer
0530   double scEtBS = ptFast(new_sc.energy(), new_sc.position(), beamSpot_->position());
0531 
0532   if (scEtBS > threshSuperClusterEt_) {
0533     switch (seed.ptr()->layer()) {
0534       case PFLayer::ECAL_BARREL:
0535         if (isOOTCollection_) {
0536           DetId seedId = new_sc.seed()->seed();
0537           EcalRecHitCollection::const_iterator seedRecHit = barrelRecHits_->find(seedId);
0538           if (!seedRecHit->checkFlag(EcalRecHit::kOutOfTime))
0539             break;
0540         }
0541         superClustersEB_->push_back(new_sc);
0542         break;
0543       case PFLayer::HGCAL:
0544       case PFLayer::ECAL_ENDCAP:
0545         if (isOOTCollection_) {
0546           DetId seedId = new_sc.seed()->seed();
0547           EcalRecHitCollection::const_iterator seedRecHit = endcapRecHits_->find(seedId);
0548           if (!seedRecHit->checkFlag(EcalRecHit::kOutOfTime))
0549             break;
0550         }
0551         superClustersEE_->push_back(new_sc);
0552         break;
0553       default:
0554         break;
0555     }
0556   }
0557 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!