File indexing completed on 2023-10-25 09:59:34
0001
0002 #include <Math/VectorUtil.h>
0003
0004 #include "RecoHI/HiEgammaAlgos/interface/EcalClusterIsoCalculator.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006
0007 #include "DataFormats/Common/interface/Handle.h"
0008 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
0009 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0010 #include "DataFormats/Math/interface/deltaR.h"
0011 #include "DataFormats/Math/interface/Vector3D.h"
0012
0013 using namespace edm;
0014 using namespace reco;
0015 using namespace std;
0016
0017 EcalClusterIsoCalculator::EcalClusterIsoCalculator(const edm::Handle<BasicClusterCollection> pEBclusters,
0018 const edm::Handle<BasicClusterCollection> pEEclusters) {
0019 if (pEBclusters.isValid())
0020 fEBclusters_ = pEBclusters.product();
0021 else
0022 fEBclusters_ = nullptr;
0023
0024 if (pEEclusters.isValid())
0025 fEEclusters_ = pEEclusters.product();
0026 else
0027 fEEclusters_ = nullptr;
0028 }
0029
0030 double EcalClusterIsoCalculator::getEcalClusterIso(const reco::SuperClusterRef cluster,
0031 const double x,
0032 const double threshold) {
0033 if (!fEBclusters_) {
0034 return -100;
0035 }
0036
0037 if (!fEEclusters_) {
0038 return -100;
0039 }
0040
0041 math::XYZVector SClusPoint(cluster->position().x(), cluster->position().y(), cluster->position().z());
0042
0043 double TotalEt = 0;
0044
0045 TotalEt = -cluster->rawEnergy() / cosh(cluster->eta());
0046
0047
0048 for (BasicClusterCollection::const_iterator iclu = fEBclusters_->begin(); iclu != fEBclusters_->end(); ++iclu) {
0049 const BasicCluster *clu = &(*iclu);
0050 math::XYZVector ClusPoint(clu->x(), clu->y(), clu->z());
0051 double eta = ClusPoint.eta();
0052
0053 double dR2 = reco::deltaR2(*clu, *cluster);
0054
0055 if (dR2 < (x * x * 0.01)) {
0056 double et = clu->energy() / cosh(eta);
0057 if (et < threshold)
0058 et = 0;
0059 TotalEt += et;
0060 }
0061 }
0062
0063 for (BasicClusterCollection::const_iterator iclu = fEEclusters_->begin(); iclu != fEEclusters_->end(); ++iclu) {
0064 const BasicCluster *clu = &(*iclu);
0065 math::XYZVector ClusPoint(clu->x(), clu->y(), clu->z());
0066 double eta = ClusPoint.eta();
0067
0068 double dR2 = reco::deltaR2(*clu, *cluster);
0069
0070 if (dR2 < (x * x * 0.01)) {
0071 double et = clu->energy() / cosh(eta);
0072 if (et < threshold)
0073 et = 0;
0074 TotalEt += et;
0075 }
0076 }
0077
0078 return TotalEt;
0079 }
0080
0081 double EcalClusterIsoCalculator::getBkgSubEcalClusterIso(const reco::SuperClusterRef cluster,
0082 const double x,
0083 double const threshold) {
0084 if (!fEBclusters_) {
0085 return -100;
0086 }
0087
0088 if (!fEEclusters_) {
0089 return -100;
0090 }
0091
0092 double SClusterEta = cluster->eta();
0093 double TotalEt = 0;
0094
0095 TotalEt = -cluster->rawEnergy() / cosh(cluster->eta());
0096
0097 for (BasicClusterCollection::const_iterator iclu = fEBclusters_->begin(); iclu != fEBclusters_->end(); ++iclu) {
0098 const BasicCluster *clu = &(*iclu);
0099 math::XYZVector ClusPoint(clu->x(), clu->y(), clu->z());
0100 double eta = ClusPoint.eta();
0101
0102 double dEta = fabs(eta - SClusterEta);
0103
0104 if (dEta < x * 0.1) {
0105 double et = clu->energy() / cosh(eta);
0106 if (et < threshold)
0107 et = 0;
0108 TotalEt += et;
0109 }
0110 }
0111
0112 for (BasicClusterCollection::const_iterator iclu = fEEclusters_->begin(); iclu != fEEclusters_->end(); ++iclu) {
0113 const BasicCluster *clu = &(*iclu);
0114 math::XYZVector ClusPoint(clu->x(), clu->y(), clu->z());
0115 double eta = ClusPoint.eta();
0116 double dEta = fabs(eta - SClusterEta);
0117
0118 if (dEta < x * 0.1) {
0119 double et = clu->energy() / cosh(eta);
0120 if (et < threshold)
0121 et = 0;
0122 TotalEt += et;
0123 }
0124 }
0125
0126 double Cx = getEcalClusterIso(cluster, x, threshold);
0127 double CCx = (Cx - TotalEt / 40.0 * x) * (1 / (1 - x / 40.));
0128
0129 return CCx;
0130 }