Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 12:50:02

0001 #ifndef CommonTools_RecoAlgos_RecoTrackSelectorBase_h
0002 #define CommonTools_RecoAlgos_RecoTrackSelectorBase_h
0003 
0004 #include "FWCore/Framework/interface/Event.h"
0005 #include "FWCore/Framework/interface/ConsumesCollector.h"
0006 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0007 #include "FWCore/Utilities/interface/InputTag.h"
0008 
0009 #include "DataFormats/TrackReco/interface/Track.h"
0010 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0011 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0012 #include "DataFormats/VertexReco/interface/Vertex.h"
0013 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0014 #include "DataFormats/Math/interface/deltaPhi.h"
0015 
0016 class RecoTrackSelectorBase {
0017 public:
0018   RecoTrackSelectorBase() {}
0019   RecoTrackSelectorBase(const edm::ParameterSet& cfg)
0020       : ptMin_(cfg.getParameter<double>("ptMin")),
0021         minRapidity_(cfg.getParameter<double>("minRapidity")),
0022         maxRapidity_(cfg.getParameter<double>("maxRapidity")),
0023         meanPhi_((cfg.getParameter<double>("minPhi") + cfg.getParameter<double>("maxPhi")) / 2.),
0024         rangePhi_((cfg.getParameter<double>("maxPhi") - cfg.getParameter<double>("minPhi")) / 2.),
0025         tip_(cfg.getParameter<double>("tip")),
0026         lip_(cfg.getParameter<double>("lip")),
0027         maxChi2_(cfg.getParameter<double>("maxChi2")),
0028         minHit_(cfg.getParameter<int>("minHit")),
0029         minPixelHit_(cfg.getParameter<int>("minPixelHit")),
0030         minLayer_(cfg.getParameter<int>("minLayer")),
0031         min3DLayer_(cfg.getParameter<int>("min3DLayer")),
0032         usePV_(false),
0033         invertRapidityCut_(cfg.getParameter<bool>("invertRapidityCut")) {
0034     const auto minPhi = cfg.getParameter<double>("minPhi");
0035     const auto maxPhi = cfg.getParameter<double>("maxPhi");
0036     if (minPhi >= maxPhi) {
0037       throw cms::Exception("Configuration")
0038           << "RecoTrackSelectorPhase: minPhi (" << minPhi << ") must be smaller than maxPhi (" << maxPhi
0039           << "). The range is constructed from minPhi to maxPhi around their average.";
0040     }
0041     if (minPhi >= M_PI) {
0042       throw cms::Exception("Configuration")
0043           << "RecoTrackSelectorPhase: minPhi (" << minPhi
0044           << ") must be smaller than PI. The range is constructed from minPhi to maxPhi around their average.";
0045     }
0046     if (maxPhi <= -M_PI) {
0047       throw cms::Exception("Configuration")
0048           << "RecoTrackSelectorPhase: maxPhi (" << maxPhi
0049           << ") must be larger than -PI. The range is constructed from minPhi to maxPhi around their average.";
0050     }
0051 
0052     for (const std::string& quality : cfg.getParameter<std::vector<std::string> >("quality"))
0053       quality_.push_back(reco::TrackBase::qualityByName(quality));
0054     for (const std::string& algorithm : cfg.getParameter<std::vector<std::string> >("algorithm"))
0055       algorithm_.push_back(reco::TrackBase::algoByName(algorithm));
0056     for (const std::string& algorithm : cfg.getParameter<std::vector<std::string> >("originalAlgorithm"))
0057       originalAlgorithm_.push_back(reco::TrackBase::algoByName(algorithm));
0058     for (const std::string& algorithm : cfg.getParameter<std::vector<std::string> >("algorithmMaskContains"))
0059       algorithmMask_.push_back(reco::TrackBase::algoByName(algorithm));
0060   }
0061 
0062   RecoTrackSelectorBase(const edm::ParameterSet& cfg, edm::ConsumesCollector& iC) : RecoTrackSelectorBase(cfg) {
0063     usePV_ = cfg.getParameter<bool>("usePV");
0064     bsSrcToken_ = iC.consumes<reco::BeamSpot>(cfg.getParameter<edm::InputTag>("beamSpot"));
0065     if (usePV_)
0066       vertexToken_ = iC.consumes<reco::VertexCollection>(cfg.getParameter<edm::InputTag>("vertexTag"));
0067   }
0068 
0069   void init(const edm::Event& event, const edm::EventSetup& es) {
0070     edm::Handle<reco::BeamSpot> beamSpot;
0071     event.getByToken(bsSrcToken_, beamSpot);
0072     vertex_ = beamSpot->position();
0073     if (!usePV_)
0074       return;
0075     edm::Handle<reco::VertexCollection> hVtx;
0076     event.getByToken(vertexToken_, hVtx);
0077     if (hVtx->empty())
0078       return;
0079     vertex_ = (*hVtx)[0].position();
0080   }
0081 
0082   bool operator()(const reco::TrackRef& tref) const { return (*this)(*tref); }
0083 
0084   bool operator()(const reco::Track& t) const { return (*this)(t, vertex_); }
0085 
0086   bool operator()(const reco::Track& t, const reco::Track::Point& vertex) const {
0087     bool quality_ok = true;
0088     if (!quality_.empty()) {
0089       quality_ok = false;
0090       for (unsigned int i = 0; i < quality_.size(); ++i) {
0091         if (t.quality(quality_[i])) {
0092           quality_ok = true;
0093           break;
0094         }
0095       }
0096     }
0097 
0098     bool algo_ok = true;
0099     if (!algorithm_.empty()) {
0100       if (std::find(algorithm_.begin(), algorithm_.end(), t.algo()) == algorithm_.end())
0101         algo_ok = false;
0102     }
0103     if (!originalAlgorithm_.empty() && algo_ok) {
0104       if (std::find(originalAlgorithm_.begin(), originalAlgorithm_.end(), t.originalAlgo()) == originalAlgorithm_.end())
0105         algo_ok = false;
0106     }
0107     if (!algorithmMask_.empty() && algo_ok) {
0108       if (std::find_if(algorithmMask_.begin(),
0109                        algorithmMask_.end(),
0110                        // for some reason I have to either explicitly give the return type, or use static_cast<bool>()
0111                        [&](reco::TrackBase::TrackAlgorithm algo) -> bool { return t.algoMask()[algo]; }) ==
0112           algorithmMask_.end())
0113         algo_ok = false;
0114     }
0115 
0116     const auto dphi = deltaPhi(t.phi(), meanPhi_);
0117 
0118     auto etaOk = [&](const reco::Track& p) -> bool {
0119       float eta = p.eta();
0120       if (!invertRapidityCut_)
0121         return (eta >= minRapidity_) && (eta <= maxRapidity_);
0122       else
0123         return (eta < minRapidity_ || eta > maxRapidity_);
0124     };
0125 
0126     return ((algo_ok & quality_ok) && t.hitPattern().numberOfValidHits() >= minHit_ &&
0127             t.hitPattern().numberOfValidPixelHits() >= minPixelHit_ &&
0128             t.hitPattern().trackerLayersWithMeasurement() >= minLayer_ &&
0129             t.hitPattern().pixelLayersWithMeasurement() + t.hitPattern().numberOfValidStripLayersWithMonoAndStereo() >=
0130                 min3DLayer_ &&
0131             fabs(t.pt()) >= ptMin_ && etaOk(t) && dphi >= -rangePhi_ && dphi <= rangePhi_ &&
0132             fabs(t.dxy(vertex)) <= tip_ && fabs(t.dsz(vertex)) <= lip_ && t.normalizedChi2() <= maxChi2_);
0133   }
0134 
0135 private:
0136   double ptMin_;
0137   double minRapidity_;
0138   double maxRapidity_;
0139   double meanPhi_;
0140   double rangePhi_;
0141   double tip_;
0142   double lip_;
0143   double maxChi2_;
0144   int minHit_;
0145   int minPixelHit_;
0146   int minLayer_;
0147   int min3DLayer_;
0148   bool usePV_;
0149   bool invertRapidityCut_;
0150 
0151   edm::EDGetTokenT<reco::BeamSpot> bsSrcToken_;
0152   edm::EDGetTokenT<reco::VertexCollection> vertexToken_;
0153 
0154   std::vector<reco::TrackBase::TrackQuality> quality_;
0155   std::vector<reco::TrackBase::TrackAlgorithm> algorithm_;
0156   std::vector<reco::TrackBase::TrackAlgorithm> originalAlgorithm_;
0157   std::vector<reco::TrackBase::TrackAlgorithm> algorithmMask_;
0158 
0159   reco::Track::Point vertex_;
0160 };
0161 
0162 #endif