HiEgammaAlgos/src/TrackIsoCalculator.cc

0001 #include "RecoHI/HiEgammaAlgos/interface/TrackIsoCalculator.h"
0002
0003 #include "DataFormats/Math/interface/deltaR.h"
0004
0005 using namespace edm;
0006 using namespace reco;
0007 using namespace std;
0008
0009 TrackIsoCalculator::TrackIsoCalculator(reco::TrackCollection const& trackCollection, std::string const& trackQuality)
0010     : recCollection_{trackCollection}, trackQuality_{trackQuality} {}
0011
0012 double TrackIsoCalculator::getTrackIso(reco::Photon const& cluster,
0013                                        const double x,
0014                                        const double threshold,
0015                                        const double innerDR) {
0016   double totalPt = 0;
0017
0018   for (auto const& recTrack : recCollection_) {
0019     bool goodtrack = recTrack.quality(reco::TrackBase::qualityByName(trackQuality_));
0020     if (!goodtrack)
0021       continue;
0022
0023     double pt = recTrack.pt();
0024     double dR2 = reco::deltaR2(cluster, recTrack);
0025     if (dR2 >= (0.01 * x * x))
0026       continue;
0027     if (dR2 < innerDR * innerDR)
0028       continue;
0029     if (pt > threshold)
0030       totalPt = totalPt + pt;
0031   }
0032
0033   return totalPt;
0034 }
0035
0036 double TrackIsoCalculator::getBkgSubTrackIso(reco::Photon const& cluster,
0037                                              const double x,
0038                                              const double threshold,
0039                                              const double innerDR) {
0040   double SClusterEta = cluster.eta();
0041   double totalPt = 0.0;
0042
0043   for (auto const& recTrack : recCollection_) {
0044     bool goodtrack = recTrack.quality(reco::TrackBase::qualityByName(trackQuality_));
0045     if (!goodtrack)
0046       continue;
0047
0048     double pt = recTrack.pt();
0049     if (std::abs(recTrack.eta() - SClusterEta) >= 0.1 * x)
0050       continue;
0051     if (reco::deltaR2(cluster, recTrack) < innerDR * innerDR)
0052       continue;
0053
0054     if (pt > threshold)
0055       totalPt = totalPt + pt;
0056   }
0057
0058   double Tx = getTrackIso(cluster, x, threshold, innerDR);
0059   double CTx = (Tx - totalPt / 40.0 * x) * (1 / (1 - x / 40.));
0060
0061   return CTx;
0062 }