CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/RecoTauTag/RecoTau/src/RecoTauQualityCuts.cc

Go to the documentation of this file.
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 // Get the KF track if it exists.  Otherwise, see if it has a GSF track.
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 // Quality cut implementations
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   // For some reason, the number of hits is signed
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   //std::cout << "QCUTS LIP: track vz: " << trk->vz() <<
00072     //" diff: " << difference << " cut: " << cut << std::endl;
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 // And a set of qcuts
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 // Get the set of Q cuts for a given type (i.e. gamma)
00121 bool mapAndCutByType(const PFCandidate& cand,
00122                      const RecoTauQualityCuts::QCutFuncMap& funcMap) {
00123   // Find the cuts that for this particle type
00124   RecoTauQualityCuts::QCutFuncMap::const_iterator cuts =
00125       funcMap.find(cand.particleId());
00126   // Return false if we dont' know how to deal w/ this particle type
00127   if (cuts == funcMap.end())
00128     return false;
00129   else
00130     // Otherwise AND all the cuts
00131     return AND(cand, cuts->second);
00132 }
00133 
00134 }  // end qcuts implementation namespace
00135 
00136 RecoTauQualityCuts::RecoTauQualityCuts(const edm::ParameterSet &qcuts) {
00137   // Setup all of our predicates
00138   QCutFuncCollection chargedHadronCuts;
00139   QCutFuncCollection gammaCuts;
00140   QCutFuncCollection neutralHadronCuts;
00141 
00142   // Make sure there are no extra passed options
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   // Build all the QCuts for tracks
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   // The impact parameter functions are bound to our member PV, since they
00180   // need it to compute the discriminant value.
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   // Require tracks to contribute a minimum weight to the associated vertex.
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   // Build the QCuts for gammas
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   // Build QCuts for netural hadrons
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   // Check if there are any remaining unparsed QCuts
00226   if (passedOptionSet.size()) {
00227     std::string unParsedOptions;
00228     bool thereIsABadParameter = false;
00229     BOOST_FOREACH(const std::string& option, passedOptionSet) {
00230       // Workaround for HLT - TODO FIXME
00231       if (option == "useTracksInsteadOfPFHadrons") {
00232         // Crash if true - no one should have this option enabled.
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       // If we get to this point, there is a real unknown parameter
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   // Make sure there are at least some quality cuts
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   // Map our QCut collections to the particle Ids they are associated to.
00264   qcuts_[PFCandidate::h] = chargedHadronCuts;
00265   qcuts_[PFCandidate::gamma] = gammaCuts;
00266   qcuts_[PFCandidate::h0] = neutralHadronCuts;
00267   // We use the same qcuts for muons/electrons and charged hadrons.
00268   qcuts_[PFCandidate::e] = chargedHadronCuts;
00269   qcuts_[PFCandidate::mu] = chargedHadronCuts;
00270 
00271   // Build a final level predicate that works on any PFCand
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     // Set null
00305     leadTrack_ = reco::TrackBaseRef();
00306   }
00307 }
00308 
00309 
00310 }}  // end namespace reco::tau