EgammaClusterAlgos/src/BremRecoveryClusterAlgo.cc

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   // 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 (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   // make the superclusters from the Barrel clusters - Island
0027   makeIslandSuperClusters(islandClustersBarrel_v, eb_rdeta_, eb_rdphi_);
0028   // make the superclusters from the EndCap clusters - Island
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     // Does our highest energy cluster have high enough energy?
0041     if ((*currentSeed)->energy() * sin((*currentSeed)->position().theta()) < seedTransverseEnergyThreshold)
0042       break;
0043
0044     // if yes, make it a seed for a new SuperCluster:
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     // and add the matching clusters:
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 }