File indexing completed on 2023-03-17 11:17:18
0001 #include "RecoEcal/EgammaClusterAlgos/interface/BremRecoveryClusterAlgo.h"
0002
0003 #include "DataFormats/CaloRecHit/interface/CaloCluster.h"
0004
0005 reco::SuperClusterCollection BremRecoveryClusterAlgo::makeSuperClusters(reco::CaloClusterPtrVector &clustersCollection) {
0006 const float etaBorder = 1.479;
0007
0008 superclusters_v.clear();
0009
0010
0011 reco::CaloClusterPtrVector islandClustersBarrel_v;
0012 reco::CaloClusterPtrVector islandClustersEndCap_v;
0013
0014
0015 for (reco::CaloCluster_iterator it = clustersCollection.begin(); it != clustersCollection.end(); it++) {
0016 reco::CaloClusterPtr cluster_p = *it;
0017 if (cluster_p->algo() == reco::CaloCluster::island) {
0018 if (fabs(cluster_p->position().eta()) < etaBorder) {
0019 islandClustersBarrel_v.push_back(cluster_p);
0020 } else {
0021 islandClustersEndCap_v.push_back(cluster_p);
0022 }
0023 }
0024 }
0025
0026
0027 makeIslandSuperClusters(islandClustersBarrel_v, eb_rdeta_, eb_rdphi_);
0028
0029 makeIslandSuperClusters(islandClustersEndCap_v, ec_rdeta_, ec_rdphi_);
0030
0031 return superclusters_v;
0032 }
0033
0034 #include "DataFormats/Math/interface/Vector3D.h"
0035
0036 void BremRecoveryClusterAlgo::makeIslandSuperClusters(reco::CaloClusterPtrVector &clusters_v,
0037 double etaRoad,
0038 double phiRoad) {
0039 for (reco::CaloCluster_iterator currentSeed = clusters_v.begin(); currentSeed != clusters_v.end(); ++currentSeed) {
0040
0041 if ((*currentSeed)->energy() * sin((*currentSeed)->position().theta()) < seedTransverseEnergyThreshold)
0042 break;
0043
0044
0045 double energy_ = (*currentSeed)->energy();
0046 math::XYZVector position_(
0047 (*currentSeed)->position().X(), (*currentSeed)->position().Y(), (*currentSeed)->position().Z());
0048 position_ *= energy_;
0049
0050 if (verbosity <= pINFO) {
0051 std::cout << "*****************************" << std::endl;
0052 std::cout << "******NEW SUPERCLUSTER*******" << std::endl;
0053 std::cout << "Seed R = " << (*currentSeed)->position().Rho() << std::endl;
0054 }
0055
0056
0057 reco::CaloClusterPtrVector constituentClusters;
0058 constituentClusters.push_back(*currentSeed);
0059 reco::CaloCluster_iterator currentCluster = currentSeed + 1;
0060 while (currentCluster != clusters_v.end()) {
0061 if (match(*currentSeed, *currentCluster, etaRoad, phiRoad)) {
0062 constituentClusters.push_back(*currentCluster);
0063 energy_ += (*currentCluster)->energy();
0064 position_ += (*currentCluster)->energy() * math::XYZVector((*currentCluster)->position().X(),
0065 (*currentCluster)->position().Y(),
0066 (*currentCluster)->position().Z());
0067 if (verbosity <= pINFO) {
0068 std::cout << "Cluster R = " << (*currentCluster)->position().Rho() << std::endl;
0069 }
0070 }
0071 ++currentCluster;
0072 }
0073
0074 position_ /= energy_;
0075
0076 if (verbosity <= pINFO) {
0077 std::cout << "Final SuperCluster R = " << position_.Rho() << std::endl;
0078 }
0079
0080 reco::SuperCluster newSuperCluster(
0081 energy_, math::XYZPoint(position_.X(), position_.Y(), position_.Z()), (*currentSeed), constituentClusters);
0082
0083 superclusters_v.push_back(newSuperCluster);
0084
0085 if (verbosity <= pINFO) {
0086 std::cout << "created a new supercluster of: " << std::endl;
0087 std::cout << "Energy = " << newSuperCluster.energy() << std::endl;
0088 std::cout << "Position in (R, phi, theta) = (" << newSuperCluster.position().Rho() << ", "
0089 << newSuperCluster.position().phi() << ", " << newSuperCluster.position().theta() << ")" << std::endl;
0090 }
0091 }
0092 clusters_v.clear();
0093 }
0094
0095 bool BremRecoveryClusterAlgo::match(reco::CaloClusterPtr seed_p,
0096 reco::CaloClusterPtr cluster_p,
0097 double dEtaMax,
0098 double dPhiMax) {
0099 math::XYZPoint clusterPosition = cluster_p->position();
0100 math::XYZPoint seedPosition = seed_p->position();
0101
0102 double dPhi = acos(cos(seedPosition.phi() - clusterPosition.phi()));
0103
0104 double dEta = fabs(seedPosition.eta() - clusterPosition.eta());
0105
0106 if (dEta > dEtaMax)
0107 return false;
0108 if (dPhi > dPhiMax)
0109 return false;
0110
0111 return true;
0112 }