HiEgammaAlgos/src/HiBremRecoveryClusterAlgo.cc

0001 #include "RecoHI/HiEgammaAlgos/interface/HiBremRecoveryClusterAlgo.h"
0002 #include "DataFormats/Math/interface/Vector3D.h"
0003 #include "DataFormats/CaloRecHit/interface/CaloCluster.h"
0004
0005 reco::SuperClusterCollection HiBremRecoveryClusterAlgo::makeSuperClusters(
0006     reco::CaloClusterPtrVector &clustersCollection) {
0007   const float etaBorder = 1.479;
0008
0009   superclusters_v.clear();
0010
0011   // create vectors of references to clusters of a specific origin...
0012   reco::CaloClusterPtrVector islandClustersBarrel_v;
0013   reco::CaloClusterPtrVector islandClustersEndCap_v;
0014
0015   // ...and populate them:
0016   for (reco::CaloCluster_iterator it = clustersCollection.begin(); it != clustersCollection.end(); it++) {
0017     reco::CaloClusterPtr cluster_p = *it;
0018     if (cluster_p->algo() == reco::CaloCluster::island) {
0019       if (verbosity <= pINFO) {
0020         std::cout << "Basic Cluster: (eta,phi,energy) = " << cluster_p->position().eta() << " "
0021                   << cluster_p->position().phi() << " " << cluster_p->energy() << std::endl;
0022       }
0023
0024       // if the basic cluster pass the energy threshold -> fill it into the list
0025       if (fabs(cluster_p->position().eta()) < etaBorder) {
0026         if (cluster_p->energy() > BarrelBremEnergyThreshold)
0027           islandClustersBarrel_v.push_back(cluster_p);
0028       } else {
0029         if (cluster_p->energy() > EndcapBremEnergyThreshold)
0030           islandClustersEndCap_v.push_back(cluster_p);
0031       }
0032     }
0033   }
0034
0035   // make the superclusters from the Barrel clusters - Island
0036   makeIslandSuperClusters(islandClustersBarrel_v, eb_rdeta_, eb_rdphi_);
0037
0038   // make the superclusters from the EndCap clusters - Island
0039   makeIslandSuperClusters(islandClustersEndCap_v, ec_rdeta_, ec_rdphi_);
0040
0041   return superclusters_v;
0042 }
0043
0044 void HiBremRecoveryClusterAlgo::makeIslandSuperClusters(reco::CaloClusterPtrVector &clusters_v,
0045                                                         double etaRoad,
0046                                                         double phiRoad) {
0047   // Vector of usedSeedEnergy, use the seed energy to check if this cluster is used.
0048   std::vector<double> usedSeedEnergy;
0049   usedSeedEnergy.clear();
0050
0051   // Main brem recovery loop
0052   for (reco::CaloCluster_iterator currentSeed = clusters_v.begin(); currentSeed != clusters_v.end(); ++currentSeed) {
0053     if (verbosity <= pINFO) {
0054       std::cout << "Current Cluster: " << (*currentSeed)->energy() << " "
0055                 << (std::find(usedSeedEnergy.begin(), usedSeedEnergy.end(), (*currentSeed)->energy()) !=
0056                     usedSeedEnergy.end())
0057                 << std::endl;
0058     }
0059
0060     // check this seed was not already used
0061     if (std::find(usedSeedEnergy.begin(), usedSeedEnergy.end(), (*currentSeed)->energy()) != usedSeedEnergy.end())
0062       continue;
0063
0064     // Does our highest energy cluster have high enough energy? If not, continue instead of break to be robust
0065     if ((*currentSeed)->energy() * sin((*currentSeed)->position().theta()) < seedTransverseEnergyThreshold)
0066       continue;
0067
0068     // if yes, make it a seed for a new SuperCluster, the position of the SC is calculated by energy weighted sum:
0069     double energy_ = (*currentSeed)->energy();
0070     math::XYZVector position_(
0071         (*currentSeed)->position().X(), (*currentSeed)->position().Y(), (*currentSeed)->position().Z());
0072     position_ *= energy_;
0073     usedSeedEnergy.push_back((*currentSeed)->energy());
0074
0075     // Printout if verbose
0076     if (verbosity <= pINFO) {
0077       std::cout << "*****************************" << std::endl;
0078       std::cout << "******NEW SUPERCLUSTER*******" << std::endl;
0079       std::cout << "Seed R = " << (*currentSeed)->position().Rho() << std::endl;
0080       std::cout << "Seed Et = " << (*currentSeed)->energy() * sin((*currentSeed)->position().theta()) << std::endl;
0081     }
0082
0083     // and add the matching clusters:
0084     reco::CaloClusterPtrVector constituentClusters;
0085     constituentClusters.push_back(*currentSeed);
0086     reco::CaloCluster_iterator currentCluster = currentSeed + 1;
0087
0088     // loop over the clusters
0089     while (currentCluster != clusters_v.end()) {
0090       // Print out the basic clusters
0091       if (verbosity <= pINFO) {
0092         std::cout << "->Cluster: " << (*currentCluster)->energy() << " Used = "
0093                   << (std::find(usedSeedEnergy.begin(), usedSeedEnergy.end(), (*currentCluster)->energy()) !=
0094                       usedSeedEnergy.end())
0095                   << " Matched = " << match(*currentSeed, *currentCluster, etaRoad, phiRoad) << std::endl;
0096       }
0097
0098       // if it's in the search window, and unused
0099       if (match(*currentSeed, *currentCluster, etaRoad, phiRoad) &&
0100           (std::find(usedSeedEnergy.begin(), usedSeedEnergy.end(), (*currentCluster)->energy()) ==
0101            usedSeedEnergy.end())) {
0102         // Add basic cluster
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         // Add the cluster to the used list
0109         usedSeedEnergy.push_back((*currentCluster)->energy());
0110
0111         if (verbosity <= pINFO) {
0112           std::cout << "Cluster R = " << (*currentCluster)->position().Rho() << std::endl;
0113         }
0114       }
0115       ++currentCluster;
0116     }
0117
0118     // Calculate the final position
0119     position_ /= energy_;
0120
0121     if (verbosity <= pINFO) {
0122       std::cout << "Final SuperCluster R = " << position_.Rho() << std::endl;
0123     }
0124
0125     // Add the supercluster to the new collection
0126     reco::SuperCluster newSuperCluster(
0127         energy_, math::XYZPoint(position_.X(), position_.Y(), position_.Z()), (*currentSeed), constituentClusters);
0128
0129     superclusters_v.push_back(newSuperCluster);
0130
0131     if (verbosity <= pINFO) {
0132       std::cout << "created a new supercluster of: " << std::endl;
0133       std::cout << "Energy = " << newSuperCluster.energy() << std::endl;
0134       std::cout << "Position in (R, phi, theta, eta) = (" << newSuperCluster.position().Rho() << ", "
0135                 << newSuperCluster.position().phi() << ", " << newSuperCluster.position().theta() << ", "
0136                 << newSuperCluster.position().eta() << ")" << std::endl;
0137     }
0138   }
0139   clusters_v.clear();
0140   usedSeedEnergy.clear();
0141 }
0142
0143 bool HiBremRecoveryClusterAlgo::match(reco::CaloClusterPtr seed_p,
0144                                       reco::CaloClusterPtr cluster_p,
0145                                       double dEtaMax,
0146                                       double dPhiMax) {
0147   math::XYZPoint clusterPosition = cluster_p->position();
0148   math::XYZPoint seedPosition = seed_p->position();
0149
0150   double dPhi = acos(cos(seedPosition.phi() - clusterPosition.phi()));
0151
0152   double dEta = fabs(seedPosition.eta() - clusterPosition.eta());
0153   if (verbosity <= pINFO) {
0154     std::cout << "seed phi: " << seedPosition.phi() << " cluster phi: " << clusterPosition.phi() << " dphi = " << dPhi
0155               << " dphiMax = " << dPhiMax << std::endl;
0156     std::cout << "seed eta: " << seedPosition.eta() << " cluster eta: " << clusterPosition.eta() << " deta = " << dEta
0157               << " detaMax = " << dEtaMax << std::endl;
0158   }
0159   if (dEta > dEtaMax)
0160     return false;
0161   if (dPhi > dPhiMax)
0162     return false;
0163
0164   return true;
0165 }