HiEgammaAlgos/src/EcalClusterIsoCalculator.cc

0001 // ROOT includes
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   // Loop over barrel basic clusters
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 }