File indexing completed on 2023-03-17 11:17:19
0001 #include "RecoEcal/EgammaClusterAlgos/interface/Multi5x5BremRecoveryClusterAlgo.h"
0002 #include "RecoEcal/EgammaCoreTools/interface/BremRecoveryPhiRoadAlgo.h"
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004
0005 reco::SuperClusterCollection Multi5x5BremRecoveryClusterAlgo::makeSuperClusters(
0006 const reco::CaloClusterPtrVector& clustersCollection) {
0007 const float etaBorder = 1.479;
0008 superclusters_v.clear();
0009
0010
0011 reco::CaloClusterPtrVector islandClustersBarrel_v;
0012 reco::CaloClusterPtrVector islandClustersEndCap_v;
0013
0014
0015 for (auto const& cluster_p : clustersCollection) {
0016 if (cluster_p->algo() == reco::CaloCluster::multi5x5) {
0017 if (fabs(cluster_p->position().eta()) < etaBorder) {
0018 islandClustersBarrel_v.push_back(cluster_p);
0019 } else {
0020 islandClustersEndCap_v.push_back(cluster_p);
0021 }
0022 }
0023 }
0024
0025
0026 makeIslandSuperClusters(islandClustersBarrel_v, eb_rdeta_, eb_rdphi_);
0027
0028 makeIslandSuperClusters(islandClustersEndCap_v, ec_rdeta_, ec_rdphi_);
0029
0030 return superclusters_v;
0031 }
0032
0033 #include "DataFormats/Math/interface/Vector3D.h"
0034
0035 void Multi5x5BremRecoveryClusterAlgo::makeIslandSuperClusters(reco::CaloClusterPtrVector& clusters_v,
0036 double etaRoad,
0037 double phiRoad) {
0038 if (clusters_v.empty())
0039 return;
0040
0041 const auto clustersSize = clusters_v.size();
0042 assert(clustersSize > 0);
0043
0044 bool usedSeed[clustersSize];
0045 for (auto ic = 0U; ic < clustersSize; ++ic)
0046 usedSeed[ic] = false;
0047
0048 float eta[clustersSize], phi[clustersSize], et[clustersSize];
0049 for (auto ic = 0U; ic < clustersSize; ++ic) {
0050 eta[ic] = clusters_v[ic]->eta();
0051 phi[ic] = clusters_v[ic]->phi();
0052 et[ic] = clusters_v[ic]->energy() * sin(clusters_v[ic]->position().theta());
0053 }
0054
0055 for (auto is = 0U; is < clustersSize; ++is) {
0056
0057 if (usedSeed[is])
0058 continue;
0059 auto const& currentSeed = clusters_v[is];
0060
0061
0062
0063 if (et[is] < seedTransverseEnergyThreshold)
0064 continue;
0065
0066
0067 double energy = (currentSeed)->energy();
0068 math::XYZVector position_(
0069 (currentSeed)->position().X(), (currentSeed)->position().Y(), (currentSeed)->position().Z());
0070 position_ *= energy;
0071 usedSeed[is] = true;
0072
0073 LogTrace("EcalClusters") << "*****************************";
0074 LogTrace("EcalClusters") << "******NEW SUPERCLUSTER*******";
0075 LogTrace("EcalClusters") << "Seed R = " << (currentSeed)->position().Rho();
0076
0077 reco::CaloClusterPtrVector constituentClusters;
0078 constituentClusters.push_back(currentSeed);
0079 auto ic = is + 1;
0080
0081 while (ic < clustersSize) {
0082 auto const& currentCluster = clusters_v[ic];
0083
0084
0085
0086 if (dynamicPhiRoad_)
0087 phiRoad = phiRoadAlgo_->endcapPhiRoad(energy + (currentCluster)->energy());
0088
0089 auto dphi = [](float p1, float p2) {
0090 auto dp = std::abs(p1 - p2);
0091 if (dp > float(M_PI))
0092 dp -= float(2 * M_PI);
0093 return std::abs(dp);
0094 };
0095
0096 auto match = [&](int i, int j) {
0097 return (dphi(phi[i], phi[j]) < phiRoad) && (std::abs(eta[i] - eta[j]) < etaRoad);
0098 };
0099
0100
0101 if (!usedSeed[ic] && match(is, ic)) {
0102
0103 constituentClusters.push_back(currentCluster);
0104 energy += (currentCluster)->energy();
0105 position_ += (currentCluster)->energy() * math::XYZVector((currentCluster)->position().X(),
0106 (currentCluster)->position().Y(),
0107 (currentCluster)->position().Z());
0108
0109 usedSeed[ic] = true;
0110 LogTrace("EcalClusters") << "Cluster R = " << (currentCluster)->position().Rho();
0111 }
0112 ++ic;
0113 }
0114
0115 position_ /= energy;
0116
0117 LogTrace("EcalClusters") << "Final SuperCluster R = " << position_.Rho();
0118
0119 reco::SuperCluster newSuperCluster(
0120 energy, math::XYZPoint(position_.X(), position_.Y(), position_.Z()), currentSeed, constituentClusters);
0121
0122 superclusters_v.push_back(newSuperCluster);
0123 LogTrace("EcalClusters") << "created a new supercluster of: ";
0124 LogTrace("EcalClusters") << "Energy = " << newSuperCluster.energy();
0125 LogTrace("EcalClusters") << "Position in (R, phi, theta) = (" << newSuperCluster.position().Rho() << ", "
0126 << newSuperCluster.position().phi() << ", " << newSuperCluster.position().theta() << ")";
0127 }
0128
0129 clusters_v.clear();
0130 }