00001 #include "RecoTauTag/RecoTau/interface/RecoTauQualityCuts.h"
00002 #include "DataFormats/VertexReco/interface/Vertex.h"
00003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00004 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
00005 #include "DataFormats/TrackReco/interface/Track.h"
00006
00007 #include <boost/bind.hpp>
00008
00009 namespace reco { namespace tau {
00010
00011 namespace {
00012
00013 const reco::TrackBaseRef getTrack(const PFCandidate& cand) {
00014 if (cand.trackRef().isNonnull())
00015 return reco::TrackBaseRef(cand.trackRef());
00016 else if (cand.gsfTrackRef().isNonnull()) {
00017 return reco::TrackBaseRef(cand.gsfTrackRef());
00018 }
00019 return reco::TrackBaseRef();
00020 }
00021 }
00022
00023
00024 namespace qcuts {
00025
00026 bool ptMin(const PFCandidate& cand, double cut) {
00027 return cand.pt() > cut;
00028 }
00029
00030 bool etMin(const PFCandidate& cand, double cut) {
00031 return cand.et() > cut;
00032 }
00033
00034 bool trkPixelHits(const PFCandidate& cand, int cut) {
00035
00036 TrackBaseRef trk = getTrack(cand);
00037 if (!trk) return false;
00038 return trk->hitPattern().numberOfValidPixelHits() >= cut;
00039 }
00040
00041 bool trkTrackerHits(const PFCandidate& cand, int cut) {
00042 TrackBaseRef trk = getTrack(cand);
00043 if (!trk) return false;
00044 return trk->hitPattern().numberOfValidHits() >= cut;
00045 }
00046
00047 bool trkTransverseImpactParameter(const PFCandidate& cand,
00048 const reco::VertexRef* pv,
00049 double cut) {
00050 if (pv->isNull()) {
00051 edm::LogError("QCutsNoPrimaryVertex") << "Primary vertex Ref in " <<
00052 "RecoTauQualityCuts is invalid. - trkTransverseImpactParameter";
00053 return false;
00054 }
00055 TrackBaseRef trk = getTrack(cand);
00056 if (!trk) return false;
00057 return std::abs(trk->dxy((*pv)->position())) <= cut;
00058 }
00059
00060 bool trkLongitudinalImpactParameter(const PFCandidate& cand,
00061 const reco::VertexRef* pv,
00062 double cut) {
00063 if (pv->isNull()) {
00064 edm::LogError("QCutsNoPrimaryVertex") << "Primary vertex Ref in " <<
00065 "RecoTauQualityCuts is invalid. - trkLongitudinalImpactParameter";
00066 return false;
00067 }
00068 TrackBaseRef trk = getTrack(cand);
00069 if (!trk) return false;
00070 double difference = std::abs(trk->dz((*pv)->position()));
00071
00072
00073 return difference <= cut;
00074 }
00075
00077 bool trkLongitudinalImpactParameterWrtTrack(const PFCandidate& cand,
00078 const reco::TrackBaseRef* trk, double cut) {
00079 if (trk->isNull()) {
00080 edm::LogError("QCutsNoValidLeadTrack") << "Lead track Ref in " <<
00081 "RecoTauQualityCuts is invalid. - trkLongitudinalImpactParameterWrtTrack";
00082 return false;
00083 }
00084 TrackBaseRef candTrk = getTrack(cand);
00085 if (!candTrk) return false;
00086 double difference = std::abs((*trk)->vz() - candTrk->vz());
00087 return difference <= cut;
00088 }
00089
00090
00091 bool minTrackVertexWeight(const PFCandidate& cand, const reco::VertexRef* pv,
00092 double cut) {
00093 if (pv->isNull()) {
00094 edm::LogError("QCutsNoPrimaryVertex") << "Primary vertex Ref in " <<
00095 "RecoTauQualityCuts is invalid. - minTrackVertexWeight";
00096 return false;
00097 }
00098 TrackBaseRef trk = getTrack(cand);
00099 if (!trk) return false;
00100 double weight = (*pv)->trackWeight(trk);
00101 return weight >= cut;
00102 }
00103
00104 bool trkChi2(const PFCandidate& cand, double cut) {
00105 TrackBaseRef trk = getTrack(cand);
00106 if (!trk) return false;
00107 return trk->normalizedChi2() <= cut;
00108 }
00109
00110
00111 bool AND(const PFCandidate& cand,
00112 const RecoTauQualityCuts::QCutFuncCollection& cuts) {
00113 BOOST_FOREACH(const RecoTauQualityCuts::QCutFunc& func, cuts) {
00114 if (!func(cand))
00115 return false;
00116 }
00117 return true;
00118 }
00119
00120
00121 bool mapAndCutByType(const PFCandidate& cand,
00122 const RecoTauQualityCuts::QCutFuncMap& funcMap) {
00123
00124 RecoTauQualityCuts::QCutFuncMap::const_iterator cuts =
00125 funcMap.find(cand.particleId());
00126
00127 if (cuts == funcMap.end())
00128 return false;
00129 else
00130
00131 return AND(cand, cuts->second);
00132 }
00133
00134 }
00135
00136 RecoTauQualityCuts::RecoTauQualityCuts(const edm::ParameterSet &qcuts) {
00137
00138 QCutFuncCollection chargedHadronCuts;
00139 QCutFuncCollection gammaCuts;
00140 QCutFuncCollection neutralHadronCuts;
00141
00142
00143 std::set<std::string> passedOptionSet;
00144 std::vector<std::string> passedOptions = qcuts.getParameterNames();
00145
00146 BOOST_FOREACH(const std::string& option, passedOptions) {
00147 passedOptionSet.insert(option);
00148 }
00149
00150
00151 if (qcuts.exists("minTrackPt")) {
00152 chargedHadronCuts.push_back(
00153 boost::bind(qcuts::ptMin, _1,
00154 qcuts.getParameter<double>("minTrackPt")));
00155 passedOptionSet.erase("minTrackPt");
00156 }
00157
00158 if (qcuts.exists("maxTrackChi2")) {
00159 chargedHadronCuts.push_back(
00160 boost::bind(qcuts::trkChi2, _1,
00161 qcuts.getParameter<double>("maxTrackChi2")));
00162 passedOptionSet.erase("maxTrackChi2");
00163 }
00164
00165 if (qcuts.exists("minTrackPixelHits")) {
00166 chargedHadronCuts.push_back(boost::bind(
00167 qcuts::trkPixelHits, _1,
00168 qcuts.getParameter<uint32_t>("minTrackPixelHits")));
00169 passedOptionSet.erase("minTrackPixelHits");
00170 }
00171
00172 if (qcuts.exists("minTrackHits")) {
00173 chargedHadronCuts.push_back(boost::bind(
00174 qcuts::trkTrackerHits, _1,
00175 qcuts.getParameter<uint32_t>("minTrackHits")));
00176 passedOptionSet.erase("minTrackHits");
00177 }
00178
00179
00180
00181 if (qcuts.exists("maxTransverseImpactParameter")) {
00182 chargedHadronCuts.push_back(boost::bind(
00183 qcuts::trkTransverseImpactParameter, _1, &pv_,
00184 qcuts.getParameter<double>("maxTransverseImpactParameter")));
00185 passedOptionSet.erase("maxTransverseImpactParameter");
00186 }
00187
00188 if (qcuts.exists("maxDeltaZ")) {
00189 chargedHadronCuts.push_back(boost::bind(
00190 qcuts::trkLongitudinalImpactParameter, _1, &pv_,
00191 qcuts.getParameter<double>("maxDeltaZ")));
00192 passedOptionSet.erase("maxDeltaZ");
00193 }
00194
00195 if (qcuts.exists("maxDeltaZToLeadTrack")) {
00196 chargedHadronCuts.push_back(boost::bind(
00197 qcuts::trkLongitudinalImpactParameterWrtTrack, _1, &leadTrack_,
00198 qcuts.getParameter<double>("maxDeltaZToLeadTrack")));
00199 passedOptionSet.erase("maxDeltaZToLeadTrack");
00200 }
00201
00202
00203 if (qcuts.exists("minTrackVertexWeight")) {
00204 chargedHadronCuts.push_back(boost::bind(
00205 qcuts::minTrackVertexWeight, _1, &pv_,
00206 qcuts.getParameter<double>("minTrackVertexWeight")));
00207 passedOptionSet.erase("minTrackVertexWeight");
00208 }
00209
00210
00211 if (qcuts.exists("minGammaEt")) {
00212 gammaCuts.push_back(boost::bind(
00213 qcuts::etMin, _1, qcuts.getParameter<double>("minGammaEt")));
00214 passedOptionSet.erase("minGammaEt");
00215 }
00216
00217
00218 if (qcuts.exists("minNeutralHadronEt")) {
00219 neutralHadronCuts.push_back(boost::bind(
00220 qcuts::etMin, _1,
00221 qcuts.getParameter<double>("minNeutralHadronEt")));
00222 passedOptionSet.erase("minNeutralHadronEt");
00223 }
00224
00225
00226 if (passedOptionSet.size()) {
00227 std::string unParsedOptions;
00228 bool thereIsABadParameter = false;
00229 BOOST_FOREACH(const std::string& option, passedOptionSet) {
00230
00231 if (option == "useTracksInsteadOfPFHadrons") {
00232
00233 if (qcuts.getParameter<bool>("useTracksInsteadOfPFHadrons")) {
00234 throw cms::Exception("DontUseTracksInQcuts")
00235 << "The obsolete exception useTracksInsteadOfPFHadrons "
00236 << "is set to true in the quality cut config." << std::endl;
00237 }
00238 continue;
00239 }
00240
00241
00242 thereIsABadParameter = true;
00243
00244 unParsedOptions += option;
00245 unParsedOptions += "\n";
00246 }
00247 if (thereIsABadParameter) {
00248 throw cms::Exception("BadQualityCutConfig")
00249 << " The PSet passed to the RecoTauQualityCuts class had"
00250 << " the following unrecognized options: " << std::endl
00251 << unParsedOptions;
00252 }
00253 }
00254
00255
00256 size_t nCuts = chargedHadronCuts.size() + gammaCuts.size()
00257 + neutralHadronCuts.size();
00258 if (!nCuts) {
00259 throw cms::Exception("BadQualityCutConfig")
00260 << " No options were passed to the quality cut class!" << std::endl;
00261 }
00262
00263
00264 qcuts_[PFCandidate::h] = chargedHadronCuts;
00265 qcuts_[PFCandidate::gamma] = gammaCuts;
00266 qcuts_[PFCandidate::h0] = neutralHadronCuts;
00267
00268 qcuts_[PFCandidate::e] = chargedHadronCuts;
00269 qcuts_[PFCandidate::mu] = chargedHadronCuts;
00270
00271
00272 predicate_ = boost::bind(qcuts::mapAndCutByType, _1, boost::cref(qcuts_));
00273 }
00274
00275 std::pair<edm::ParameterSet, edm::ParameterSet> factorizePUQCuts(
00276 const edm::ParameterSet& input) {
00277
00278 edm::ParameterSet puCuts;
00279 edm::ParameterSet nonPUCuts;
00280
00281 std::vector<std::string> inputNames = input.getParameterNames();
00282 BOOST_FOREACH(const std::string& cut, inputNames) {
00283 if (cut == "minTrackVertexWeight" || cut == "maxDeltaZ"
00284 || cut == "maxDeltaZToLeadTrack") {
00285 puCuts.copyFrom(input, cut);
00286 } else {
00287 nonPUCuts.copyFrom(input, cut);
00288 }
00289 }
00290 return std::make_pair(puCuts, nonPUCuts);
00291 }
00292
00293
00294 void RecoTauQualityCuts::setLeadTrack(
00295 const reco::PFCandidate& leadCand) const {
00296 leadTrack_ = getTrack(leadCand);
00297 }
00298
00299 void RecoTauQualityCuts::setLeadTrack(
00300 const reco::PFCandidateRef& leadCand) const {
00301 if (leadCand.isNonnull()) {
00302 leadTrack_ = getTrack(*leadCand);
00303 } else {
00304
00305 leadTrack_ = reco::TrackBaseRef();
00306 }
00307 }
00308
00309
00310 }}