EgammaIsolationAlgos/src/PhotonTkIsolation.cc

0001 //*****************************************************************************
0002 // File:      PhotonTkIsolation.cc
0003 // ----------------------------------------------------------------------------
0004 // OrigAuth:  Matthias Mozer
0005 // Institute: IIHE-VUB
0006 //=============================================================================
0007 //*****************************************************************************
0008
0009 #include "RecoEgamma/EgammaIsolationAlgos/interface/PhotonTkIsolation.h"
0010 #include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaTrackSelector.h"
0011
0012 PhotonTkIsolation::PhotonTkIsolation(float extRadius,
0013                                      float intRadiusBarrel,
0014                                      float intRadiusEndcap,
0015                                      float stripBarrel,
0016                                      float stripEndcap,
0017                                      float etLow,
0018                                      float lip,
0019                                      float drb,
0020                                      const reco::TrackCollection* trackCollection,
0021                                      reco::TrackBase::Point beamPoint,
0022                                      const std::string& dzOptionString)
0023     : extRadius2_(extRadius * extRadius),
0024       intRadiusBarrel2_(intRadiusBarrel * intRadiusBarrel),
0025       intRadiusEndcap2_(intRadiusEndcap * intRadiusEndcap),
0026       stripBarrel_(stripBarrel),
0027       stripEndcap_(stripEndcap),
0028       etLow_(etLow),
0029       lip_(lip),
0030       drb_(drb),
0031       trackCollection_(trackCollection),
0032       beamPoint_(beamPoint) {
0033   setDzOption(dzOptionString);
0034 }
0035
0036 void PhotonTkIsolation::setDzOption(const std::string& s) {
0037   if (!s.compare("dz"))
0038     dzOption_ = egammaisolation::EgammaTrackSelector::dz;
0039   else if (!s.compare("vz"))
0040     dzOption_ = egammaisolation::EgammaTrackSelector::vz;
0041   else if (!s.compare("bs"))
0042     dzOption_ = egammaisolation::EgammaTrackSelector::bs;
0043   else if (!s.compare("vtx"))
0044     dzOption_ = egammaisolation::EgammaTrackSelector::vtx;
0045   else
0046     dzOption_ = egammaisolation::EgammaTrackSelector::dz;
0047 }
0048
0049 PhotonTkIsolation::~PhotonTkIsolation() {}
0050
0051 // unified acces to isolations
0052 std::pair<int, float> PhotonTkIsolation::getIso(const reco::Candidate* photon) const {
0053   int counter = 0;
0054   float ptSum = 0.;
0055
0056   //Take the photon position
0057   float photonEta = photon->eta();
0058
0059   //loop over tracks
0060   for (reco::TrackCollection::const_iterator trItr = trackCollection_->begin(); trItr != trackCollection_->end();
0061        ++trItr) {
0062     //check z-distance of vertex
0063     float dzCut = 0;
0064     switch (dzOption_) {
0065       case egammaisolation::EgammaTrackSelector::dz:
0066         dzCut = fabs((*trItr).dz() - photon->vertex().z());
0067         break;
0068       case egammaisolation::EgammaTrackSelector::vz:
0069         dzCut = fabs((*trItr).vz() - photon->vertex().z());
0070         break;
0071       case egammaisolation::EgammaTrackSelector::bs:
0072         dzCut = fabs((*trItr).dz(beamPoint_) - photon->vertex().z());
0073         break;
0074       case egammaisolation::EgammaTrackSelector::vtx:
0075         dzCut = fabs((*trItr).dz(photon->vertex()));
0076         break;
0077       default:
0078         dzCut = fabs((*trItr).vz() - photon->vertex().z());
0079         break;
0080     }
0081     if (dzCut > lip_)
0082       continue;
0083
0084     float this_pt = (*trItr).pt();
0085     if (this_pt < etLow_)
0086       continue;
0087     if (fabs((*trItr).dxy(beamPoint_)) > drb_)
0088       continue;  // only consider tracks from the main vertex
0089     float dr2 = reco::deltaR2(*trItr, *photon);
0090     float deta = (*trItr).eta() - photonEta;
0091     if (fabs(photonEta) < 1.479) {
0092       if (dr2 < extRadius2_ && dr2 >= intRadiusBarrel2_ && fabs(deta) >= stripBarrel_) {
0093         ++counter;
0094         ptSum += this_pt;
0095       }
0096     } else {
0097       if (dr2 < extRadius2_ && dr2 >= intRadiusEndcap2_ && fabs(deta) >= stripEndcap_) {
0098         ++counter;
0099         ptSum += this_pt;
0100       }
0101     }
0102
0103   }  //end loop over tracks
0104
0105   std::pair<int, float> retval;
0106   retval.first = counter;
0107   retval.second = ptSum;
0108   return retval;
0109 }