ElectronTools/src/SuperClusterHelper.cc

0001 #include "EgammaAnalysis/ElectronTools/interface/SuperClusterHelper.h"
0002 #include "Geometry/CaloTopology/interface/CaloTopology.h"
0003 #include "FWCore/Utilities/interface/isFinite.h"
0004 #include "RecoEgamma/EgammaTools/interface/EcalClusterLocal.h"
0005
0006 SuperClusterHelper::SuperClusterHelper(const reco::GsfElectron* electron,
0007                                        const EcalRecHitCollection* rechits,
0008                                        const CaloTopology* topo,
0009                                        const CaloGeometry* geom) {
0010   theElectron_ = electron;
0011   rechits_ = rechits;
0012   seedCluster_ = &(*(electron->superCluster()->seed()));
0013   theSuperCluster_ = &(*electron->superCluster());
0014
0015   eSubClusters_ = 0.;
0016   // Store subclusters
0017   reco::CaloCluster_iterator itscl = theSuperCluster_->clustersBegin();
0018   reco::CaloCluster_iterator itsclE = theSuperCluster_->clustersEnd();
0019   for (; itscl < itsclE; ++itscl) {
0020     if ((*itscl) == electron->superCluster()->seed())
0021       continue;  // skip seed cluster
0022     theBasicClusters_.push_back(&(**itscl));
0023     eSubClusters_ += (*itscl)->energy();
0024   }
0025   // sort subclusters
0026   sort(theBasicClusters_.begin(), theBasicClusters_.end(), SuperClusterHelper::sortClusters);
0027   // Add seed cluster at the beginning
0028   theBasicClusters_.insert(theBasicClusters_.begin(), seedCluster_);
0029
0030   nBasicClusters_ = theBasicClusters_.size();
0031
0032   // Store ES clusters
0033   eESClusters_ = 0.;
0034   itscl = theSuperCluster_->preshowerClustersBegin();
0035   itsclE = theSuperCluster_->preshowerClustersEnd();
0036   for (; itscl < itsclE; ++itscl) {
0037     theESClusters_.push_back(&(**itscl));
0038     eESClusters_ += (*itscl)->energy();
0039   }
0040   // sort ES clusters
0041   sort(theESClusters_.begin(), theESClusters_.end(), SuperClusterHelper::sortClusters);
0042
0043   nESClusters_ = theESClusters_.size();
0044
0045   topology_ = topo;
0046   geometry_ = geom;
0047   barrel_ = electron->isEB();
0048   covComputed_ = false;
0049   localCoordinatesComputed_ = false;
0050 }
0051
0052 SuperClusterHelper::SuperClusterHelper(const pat::Electron* electron,
0053                                        const EcalRecHitCollection* rechits,
0054                                        const CaloTopology* topo,
0055                                        const CaloGeometry* geom) {
0056   theElectron_ = (const reco::GsfElectron*)electron;
0057   rechits_ = rechits;
0058   //  for(unsigned ir=0; ir<rechits_->size();++ir) {
0059   //    std::cout << "RecHit " << (*rechits_)[ir].id().rawId() << " " << (*rechits_)[ir] << std::endl;
0060   //  }
0061   // Get the embedded objects
0062   theSuperCluster_ = &(*electron->superCluster());
0063   seedCluster_ = &(*(electron->seed()));
0064   const std::vector<reco::CaloCluster>& basicClusters(electron->basicClusters());
0065   nBasicClusters_ = basicClusters.size();
0066   eSubClusters_ = 0.;
0067   // Store subclusters
0068   for (unsigned ib = 0; ib < nBasicClusters_; ++ib) {
0069     if (fabs((basicClusters[ib].energy() - seedCluster_->energy()) / seedCluster_->energy()) < 1.e-5 &&
0070         fabs((basicClusters[ib].eta() - seedCluster_->eta()) / seedCluster_->eta()) < 1.e-5 &&
0071         fabs((basicClusters[ib].phi() - seedCluster_->phi()) / seedCluster_->phi()) < 1.e-5)
0072       continue;  // skip seed cluster
0073     theBasicClusters_.push_back(&basicClusters[ib]);
0074     eSubClusters_ += basicClusters[ib].energy();
0075   }
0076   // sort subclusters
0077   sort(theBasicClusters_.begin(), theBasicClusters_.end(), SuperClusterHelper::sortClusters);
0078   // Add seed cluster at the beginning
0079   theBasicClusters_.insert(theBasicClusters_.begin(), seedCluster_);
0080
0081   // Store ES clusters
0082   const std::vector<reco::CaloCluster>& esClusters(electron->preshowerClusters());
0083   nESClusters_ = esClusters.size();
0084   eESClusters_ = 0.;
0085   for (unsigned ib = 0; ib < nESClusters_; ++ib) {
0086     theESClusters_.push_back(&esClusters[ib]);
0087     eESClusters_ += esClusters[ib].energy();
0088   }
0089   // sort ES clusters
0090   sort(theESClusters_.begin(), theESClusters_.end(), SuperClusterHelper::sortClusters);
0091
0092   //  std::vector< std::pair<DetId, float> >::const_iterator it=seedCluster_->hitsAndFractions().begin();
0093   //  std::vector< std::pair<DetId, float> >::const_iterator itend=seedCluster_->hitsAndFractions().end();
0094   //  for( ; it!=itend ; ++it) {
0095   //    DetId id=it->first;
0096   //    std::cout << " Basic cluster " << id.rawId() << std::endl;
0097   //  }
0098   topology_ = topo;
0099   geometry_ = geom;
0100   barrel_ = electron->isEB();
0101   covComputed_ = false;
0102   localCoordinatesComputed_ = false;
0103 }
0104
0105 void SuperClusterHelper::computeLocalCovariances() {
0106   if (!covComputed_) {
0107     const auto& vCov_ = EcalClusterTools::localCovariances(*seedCluster_, rechits_, topology_, 4.7);
0108     covComputed_ = true;
0109
0110     spp_ = 0;
0111     if (edm::isFinite(vCov_[2]))
0112       spp_ = sqrt(vCov_[2]);
0113
0114     if (theElectron_->sigmaIetaIeta() * spp_ > 0) {
0115       sep_ = vCov_[1] / (theElectron_->sigmaIetaIeta() * spp_);
0116     } else if (vCov_[1] > 0) {
0117       sep_ = 1.0;
0118     } else {
0119       sep_ = -1.0;
0120     }
0121   }
0122 }
0123
0124 float SuperClusterHelper::spp() {
0125   computeLocalCovariances();
0126   return spp_;
0127 }
0128
0129 float SuperClusterHelper::sep() {
0130   computeLocalCovariances();
0131   return sep_;
0132 }
0133
0134 void SuperClusterHelper::localCoordinates() {
0135   if (localCoordinatesComputed_)
0136     return;
0137
0138   if (barrel_) {
0139     egammaTools::localEcalClusterCoordsEB(
0140         *seedCluster_, *geometry_, etaCrySeed_, phiCrySeed_, ietaSeed_, iphiSeed_, thetaTilt_, phiTilt_);
0141   } else {
0142     egammaTools::localEcalClusterCoordsEE(
0143         *seedCluster_, *geometry_, etaCrySeed_, phiCrySeed_, ietaSeed_, iphiSeed_, thetaTilt_, phiTilt_);
0144   }
0145   localCoordinatesComputed_ = true;
0146 }
0147
0148 float SuperClusterHelper::subClusterEnergy(unsigned i) const {
0149   return (nBasicClusters_ > i) ? theBasicClusters_[i]->energy() : 0.;
0150 }
0151
0152 float SuperClusterHelper::subClusterEta(unsigned i) const {
0153   return (nBasicClusters_ > i) ? theBasicClusters_[i]->eta() : 999.;
0154 }
0155
0156 float SuperClusterHelper::subClusterPhi(unsigned i) const {
0157   return (nBasicClusters_ > i) ? theBasicClusters_[i]->phi() : 999.;
0158 }
0159
0160 float SuperClusterHelper::subClusterEmax(unsigned i) const {
0161   return (nBasicClusters_ > i) ? EcalClusterTools::eMax(*theBasicClusters_[i], rechits_) : 0.;
0162 }
0163
0164 float SuperClusterHelper::subClusterE3x3(unsigned i) const {
0165   return (nBasicClusters_ > i) ? EcalClusterTools::e3x3(*theBasicClusters_[i], rechits_, topology_) : 0.;
0166 }
0167
0168 float SuperClusterHelper::esClusterEnergy(unsigned i) const {
0169   return (nESClusters_ > i) ? theESClusters_[i]->energy() : 0.;
0170 }
0171
0172 float SuperClusterHelper::esClusterEta(unsigned i) const {
0173   return (nESClusters_ > i) ? theESClusters_[i]->eta() : 999.;
0174 }
0175
0176 float SuperClusterHelper::esClusterPhi(unsigned i) const {
0177   return (nESClusters_ > i) ? theESClusters_[i]->phi() : 999.;
0178 }