EgammaCoreTools/src/EcalClusterPUCleaningTools.cc

0001 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0002 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterPUCleaningTools.h"
0003
0004 #include "Geometry/CaloTopology/interface/CaloTopology.h"
0005 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0006
0007 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0008
0009 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h"
0010
0011 #include "RecoEcal/EgammaCoreTools/interface/SuperClusterShapeAlgo.h"
0012
0013 EcalClusterPUCleaningTools::EcalClusterPUCleaningTools(edm::ConsumesCollector &cc,
0014                                                        const edm::InputTag &redEBRecHits,
0015                                                        const edm::InputTag &redEERecHits)
0016     : pEBRecHitsToken_(cc.consumes<EcalRecHitCollection>(redEBRecHits)),
0017       pEERecHitsToken_(cc.consumes<EcalRecHitCollection>(redEERecHits)),
0018       geometryToken_(cc.esConsumes()) {}
0019
0020 EcalClusterPUCleaningTools::~EcalClusterPUCleaningTools() {}
0021
0022 void EcalClusterPUCleaningTools::getEBRecHits(const edm::Event &ev) {
0023   edm::Handle<EcalRecHitCollection> pEBRecHits;
0024   ev.getByToken(pEBRecHitsToken_, pEBRecHits);
0025   ebRecHits_ = pEBRecHits.product();
0026 }
0027
0028 void EcalClusterPUCleaningTools::getEERecHits(const edm::Event &ev) {
0029   edm::Handle<EcalRecHitCollection> pEERecHits;
0030   ev.getByToken(pEERecHitsToken_, pEERecHits);
0031   eeRecHits_ = pEERecHits.product();
0032 }
0033
0034 reco::SuperCluster EcalClusterPUCleaningTools::CleanedSuperCluster(float xi,
0035                                                                    const reco::SuperCluster &scluster,
0036                                                                    const edm::Event &ev,
0037                                                                    const edm::EventSetup &es) {
0038   // get the geometry
0039   const auto &geometry = es.getData(geometryToken_);
0040
0041   // get the RecHits
0042   getEBRecHits(ev);
0043   getEERecHits(ev);
0044
0045   // seed basic cluster of initial SC: this will remain in the cleaned SC, by construction
0046   const reco::CaloClusterPtr &seed = scluster.seed();
0047
0048   // this should be replaced by the 5x5 around the seed; a good approx of how E_seed is defined
0049   float seedBCEnergy = (scluster.seed())->energy();
0050   float eSeed = 0.35;  // standard eSeed in EB ; see CMS IN-2010/008
0051
0052   double ClusterE = 0;  //Sum of cluster energies for supercluster.
0053   //Holders for position of this supercluster.
0054   double posX = 0;
0055   double posY = 0;
0056   double posZ = 0;
0057
0058   reco::CaloClusterPtrVector thissc;
0059
0060   // looping on basic clusters within the Supercluster
0061   for (reco::CaloCluster_iterator bcIt = scluster.clustersBegin(); bcIt != scluster.clustersEnd(); bcIt++) {
0062     // E_seed is an Et  selection on 5x1 dominoes (around which BC's are built), see CMS IN-2010/008
0063     // here such selection is implemented on the full BC around it
0064     if ((*bcIt)->energy() >
0065         sqrt(eSeed * eSeed + xi * xi * seedBCEnergy * seedBCEnergy / cosh((*bcIt)->eta()) / cosh((*bcIt)->eta()))) {
0066       ;
0067     }  // the sum only of the BC's that pass the Esee selection
0068
0069     // if BC passes dynamic selection, include it in the 'cleaned' supercluster
0070     thissc.push_back((*bcIt));
0071     // cumulate energy and position of the cleaned supercluster
0072     ClusterE += (*bcIt)->energy();
0073     posX += (*bcIt)->energy() * (*bcIt)->position().X();
0074     posY += (*bcIt)->energy() * (*bcIt)->position().Y();
0075     posZ += (*bcIt)->energy() * (*bcIt)->position().Z();
0076
0077   }  // loop on basic clusters of the original SC
0078
0079   posX /= ClusterE;
0080   posY /= ClusterE;
0081   posZ /= ClusterE;
0082
0083   // make temporary 'cleaned' supercluster in oder to compute the covariances
0084   double Epreshower = scluster.preshowerEnergy();
0085   double phiWidth = 0.;
0086   double etaWidth = 0.;
0087   reco::SuperCluster suCltmp(ClusterE, math::XYZPoint(posX, posY, posZ), seed, thissc, Epreshower, phiWidth, etaWidth);
0088
0089   // construct cluster shape to compute ieta and iphi covariances of the SC
0090   const CaloSubdetectorGeometry *geometry_p = nullptr;
0091   if (seed->seed().det() == DetId::Ecal && seed->seed().subdetId() == EcalBarrel) {
0092     geometry_p = geometry.getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
0093     SuperClusterShapeAlgo SCShape(ebRecHits_, geometry_p);
0094     SCShape.Calculate_Covariances(suCltmp);
0095     phiWidth = SCShape.phiWidth();
0096     etaWidth = SCShape.etaWidth();
0097   } else if (seed->seed().det() == DetId::Ecal && seed->seed().subdetId() == EcalEndcap) {
0098     geometry_p = geometry.getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
0099     SuperClusterShapeAlgo SCShape(eeRecHits_, geometry_p);
0100     SCShape.Calculate_Covariances(suCltmp);
0101     phiWidth = SCShape.phiWidth();
0102     etaWidth = SCShape.etaWidth();
0103   } else {
0104     edm::LogError("SeedError")
0105         << "The seed crystal of this SC is neither in EB nor in EE. This is a problem. Bailing out";
0106     assert(-1);
0107   }
0108
0109   // return the cleaned supercluster SCluster, with covariances updated
0110   return reco::SuperCluster(ClusterE, math::XYZPoint(posX, posY, posZ), seed, thissc, Epreshower, phiWidth, etaWidth);
0111 }