EgammaClusterAlgos/src/Multi5x5BremRecoveryClusterAlgo.cc

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   // create vectors of references to clusters of a specific origin...
0011   reco::CaloClusterPtrVector islandClustersBarrel_v;
0012   reco::CaloClusterPtrVector islandClustersEndCap_v;
0013
0014   // ...and populate them:
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   // make the superclusters from the Barrel clusters - Island
0026   makeIslandSuperClusters(islandClustersBarrel_v, eb_rdeta_, eb_rdphi_);
0027   // make the superclusters from the EndCap clusters - Island
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     // check this seed was not already used
0057     if (usedSeed[is])
0058       continue;
0059     auto const& currentSeed = clusters_v[is];
0060
0061     // Does our highest energy cluster have high enough energy?
0062     // changed this to continue from break (to be robust against the order of sorting of the seed clusters)
0063     if (et[is] < seedTransverseEnergyThreshold)
0064       continue;
0065
0066     // if yes, make it a seed for a new SuperCluster:
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       // if dynamic phi road is enabled then compute the phi road for a cluster
0085       // of energy of existing clusters + the candidate cluster.
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       // does the cluster match the phi road for this candidate supercluster
0101       if (!usedSeed[ic] && match(is, ic)) {
0102         // add basic cluster to supercluster constituents
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         // remove cluster from vector of available clusters
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 }