EgammaIsolationAlgos/src/EgammaEcalIsolation.cc

0001 //*****************************************************************************
0002 // File:      EgammaEcalIsolation.cc
0003 // ----------------------------------------------------------------------------
0004 // OrigAuth:  Gilles De Lentdecker
0005 // Institute: IIHE-ULB
0006 //=============================================================================
0007 //*****************************************************************************
0008
0009 //C++ includes
0010 #include <vector>
0011 #include <functional>
0012
0013 //ROOT includes
0014 #include <Math/VectorUtil.h>
0015
0016 //CMSSW includes
0017 #include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaEcalIsolation.h"
0018 #include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
0019
0020 using namespace ROOT::Math::VectorUtil;
0021
0022 EgammaEcalIsolation::EgammaEcalIsolation(double extRadius,
0023                                          double etLow,
0024                                          const reco::BasicClusterCollection* basicClusterCollection,
0025                                          const reco::SuperClusterCollection* superClusterCollection)
0026     : etMin(etLow),
0027       conesize(extRadius),
0028       basicClusterCollection_(basicClusterCollection),
0029       superClusterCollection_(superClusterCollection) {}
0030
0031 EgammaEcalIsolation::~EgammaEcalIsolation() {}
0032
0033 double EgammaEcalIsolation::getEcalEtSum(const reco::Candidate* candidate) {
0034   double ecalIsol = 0.;
0035   reco::SuperClusterRef sc = candidate->get<reco::SuperClusterRef>();
0036   math::XYZVector position(sc.get()->position().x(), sc.get()->position().y(), sc.get()->position().z());
0037
0038   // match the photon hybrid supercluster with those with Algo==0 (island)
0039   double delta1 = 1000.;
0040   double deltacur = 1000.;
0041   const reco::SuperCluster* matchedsupercluster = nullptr;
0042   bool MATCHEDSC = false;
0043
0044   for (reco::SuperClusterCollection::const_iterator scItr = superClusterCollection_->begin();
0045        scItr != superClusterCollection_->end();
0046        ++scItr) {
0047     const reco::SuperCluster* supercluster = &(*scItr);
0048
0049     math::XYZVector currentPosition(
0050         supercluster->position().x(), supercluster->position().y(), supercluster->position().z());
0051
0052     if (supercluster->seed()->algo() == 0) {
0053       deltacur = DeltaR(currentPosition, position);
0054
0055       if (deltacur < delta1) {
0056         delta1 = deltacur;
0057         matchedsupercluster = supercluster;
0058         MATCHEDSC = true;
0059       }
0060     }
0061   }
0062
0063   const reco::BasicCluster* cluster = nullptr;
0064
0065   //loop over basic clusters
0066   for (reco::BasicClusterCollection::const_iterator cItr = basicClusterCollection_->begin();
0067        cItr != basicClusterCollection_->end();
0068        ++cItr) {
0069     cluster = &(*cItr);
0070     //    double ebc_bcchi2 = cluster->chi2();
0071     int ebc_bcalgo = cluster->algo();
0072     double ebc_bce = cluster->energy();
0073     double ebc_bceta = cluster->eta();
0074     double ebc_bcet = ebc_bce * sin(2 * atan(exp(ebc_bceta)));
0075     double newDelta = 0.;
0076
0077     if (ebc_bcet > etMin && ebc_bcalgo == 0) {
0078       //    if (ebc_bcchi2 < 30.) {
0079
0080       if (MATCHEDSC) {
0081         bool inSuperCluster = false;
0082
0083         reco::CaloCluster_iterator theEclust = matchedsupercluster->clustersBegin();
0084         // loop over the basic clusters of the matched supercluster
0085         for (; theEclust != matchedsupercluster->clustersEnd(); theEclust++) {
0086           if ((**theEclust) == (*cluster))
0087             inSuperCluster = true;
0088         }
0089         if (!inSuperCluster) {
0090           math::XYZVector basicClusterPosition(
0091               cluster->position().x(), cluster->position().y(), cluster->position().z());
0092           newDelta = DeltaR(basicClusterPosition, position);
0093           if (newDelta < conesize) {
0094             ecalIsol += ebc_bcet;
0095           }
0096         }
0097       }
0098       //      } // matches ebc_bcchi2
0099     }  // matches ebc_bcet && ebc_bcalgo
0100   }
0101
0102   //  std::cout << "Will return ecalIsol = " << ecalIsol << std::endl;
0103   return ecalIsol;
0104 }