Line Code
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154
#include "CommonTools/RecoAlgos/interface/TrackWithVertexSelector.h"
#include "FWCore/Utilities/interface/isFinite.h"
//
// constructors and destructor
//

namespace {
  constexpr float fakeBeamSpotTimeWidth = 0.175f;
}

TrackWithVertexSelector::TrackWithVertexSelector(const edm::ParameterSet &iConfig, edm::ConsumesCollector &iC)
    : numberOfValidHits_(iConfig.getParameter<uint32_t>("numberOfValidHits")),
      numberOfValidPixelHits_(iConfig.getParameter<uint32_t>("numberOfValidPixelHits")),
      numberOfLostHits_(iConfig.getParameter<uint32_t>("numberOfLostHits")),
      normalizedChi2_(iConfig.getParameter<double>("normalizedChi2")),
      ptMin_(iConfig.getParameter<double>("ptMin")),
      ptMax_(iConfig.getParameter<double>("ptMax")),
      etaMin_(iConfig.getParameter<double>("etaMin")),
      etaMax_(iConfig.getParameter<double>("etaMax")),
      dzMax_(iConfig.getParameter<double>("dzMax")),
      d0Max_(iConfig.getParameter<double>("d0Max")),
      ptErrorCut_(iConfig.getParameter<double>("ptErrorCut")),
      quality_(iConfig.getParameter<std::string>("quality")),
      nVertices_(iConfig.getParameter<bool>("useVtx") ? iConfig.getParameter<uint32_t>("nVertices") : 0),
      vertexToken_(iC.consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("vertexTag"))),
      timesToken_(iC.consumes<edm::ValueMap<float> >(iConfig.getParameter<edm::InputTag>("timesTag"))),
      timeResosToken_(iC.consumes<edm::ValueMap<float> >(iConfig.getParameter<edm::InputTag>("timeResosTag"))),
      vtxFallback_(iConfig.getParameter<bool>("vtxFallback")),
      zetaVtx_(iConfig.getParameter<double>("zetaVtx")),
      rhoVtx_(iConfig.getParameter<double>("rhoVtx")),
      nSigmaDtVertex_(iConfig.getParameter<double>("nSigmaDtVertex")) {}

TrackWithVertexSelector::~TrackWithVertexSelector() {}

void TrackWithVertexSelector::init(const edm::Event &event) {
  edm::Handle<reco::VertexCollection> hVtx;
  event.getByToken(vertexToken_, hVtx);
  vcoll_ = nVertices_ > 0 ? hVtx.product() : nullptr;

  edm::Handle<edm::ValueMap<float> > hTimes;
  event.getByToken(timesToken_, hTimes);
  timescoll_ = hTimes.isValid() ? hTimes.product() : nullptr;

  edm::Handle<edm::ValueMap<float> > hTimeResos;
  event.getByToken(timeResosToken_, hTimeResos);
  timeresoscoll_ = hTimeResos.isValid() ? hTimeResos.product() : nullptr;
}

bool TrackWithVertexSelector::testTrack(const reco::Track &t) const {
  using std::abs;
  if ((t.numberOfValidHits() >= numberOfValidHits_) &&
      (static_cast<unsigned int>(t.hitPattern().numberOfValidPixelHits()) >= numberOfValidPixelHits_) &&
      (t.numberOfLostHits() <= numberOfLostHits_) && (t.normalizedChi2() <= normalizedChi2_) &&
      (t.ptError() / t.pt() * std::max(1., t.normalizedChi2()) <= ptErrorCut_) &&
      (t.quality(t.qualityByName(quality_))) && (t.pt() >= ptMin_) && (t.pt() <= ptMax_) && (abs(t.eta()) <= etaMax_) &&
      (abs(t.eta()) >= etaMin_) && (abs(t.dz()) <= dzMax_) && (abs(t.d0()) <= d0Max_)) {
    return true;
  }
  return false;
}

bool TrackWithVertexSelector::testTrack(const reco::TrackRef &tref) const { return testTrack(*tref); }

bool TrackWithVertexSelector::testVertices(const reco::Track &t, const reco::VertexCollection &vtxs) const {
  bool ok = false;
  if (!vtxs.empty()) {
    unsigned int tested = 1;
    for (reco::VertexCollection::const_iterator it = vtxs.begin(), ed = vtxs.end(); it != ed; ++it) {
      if ((std::abs(t.dxy(it->position())) < rhoVtx_) && (std::abs(t.dz(it->position())) < zetaVtx_)) {
        ok = true;
        break;
      }
      if (tested++ >= nVertices_)
        break;
    }
  } else if (vtxFallback_) {
    return ((std::abs(t.vertex().z()) < 15.9) && (t.vertex().Rho() < 0.2));
  }
  return ok;
}

bool TrackWithVertexSelector::testVertices(const reco::TrackRef &tref, const reco::VertexCollection &vtxs) const {
  const auto &t = *tref;
  const bool timeAvailable = timescoll_ != nullptr && timeresoscoll_ != nullptr;
  bool ok = false;
  if (!vtxs.empty()) {
    unsigned int tested = 1;
    for (reco::VertexCollection::const_iterator it = vtxs.begin(), ed = vtxs.end(); it != ed; ++it) {
      const bool useTime = timeAvailable && it->t() != 0.;
      float time = useTime ? (*timescoll_)[tref] : -1.f;
      float timeReso = useTime ? (*timeresoscoll_)[tref] : -1.f;
      timeReso = (timeReso > 1e-6 ? timeReso : fakeBeamSpotTimeWidth);

      if (edm::isNotFinite(time)) {
        time = 0.0;
        timeReso = 2.0 * fakeBeamSpotTimeWidth;
      }

      const double vtxSigmaT2 = it->tError() * it->tError();
      const double vtxTrackErr = std::sqrt(vtxSigmaT2 + timeReso * timeReso);

      if ((std::abs(t.dxy(it->position())) < rhoVtx_) && (std::abs(t.dz(it->position())) < zetaVtx_) &&
          (!useTime || (std::abs(time - it->t()) / vtxTrackErr < nSigmaDtVertex_))) {
        ok = true;
        break;
      }
      if (tested++ >= nVertices_)
        break;
    }
  } else if (vtxFallback_) {
    return ((std::abs(t.vertex().z()) < 15.9) && (t.vertex().Rho() < 0.2));
  }
  return ok;
}

bool TrackWithVertexSelector::operator()(const reco::Track &t) const {
  if (!testTrack(t))
    return false;
  if (nVertices_ == 0)
    return true;
  return testVertices(t, *vcoll_);
}

bool TrackWithVertexSelector::operator()(const reco::TrackRef &tref) const {
  if (!testTrack(tref))
    return false;
  if (nVertices_ == 0)
    return true;
  return testVertices(tref, *vcoll_);
}

void TrackWithVertexSelector::fillPSetDescription(edm::ParameterSetDescription &desc) {
  desc.add<uint32_t>("numberOfValidHits", 0);
  desc.add<uint32_t>("numberOfValidPixelHits", 0);
  desc.add<uint32_t>("numberOfLostHits", 999)->setComment("at most 999 lost hits");
  desc.add<double>("normalizedChi2", 999999.);
  desc.add<double>("ptMin", 0.3)->setComment("in GeV");
  desc.add<double>("ptMax", 500.0)->setComment("in GeV");
  desc.add<double>("etaMin", 0.0);
  desc.add<double>("etaMax", 50.0);
  desc.add<double>("dzMax", 999.)->setComment("in cm");
  desc.add<double>("d0Max", 999.)->setComment("in cm");
  desc.add<double>("ptErrorCut", 0.2)->setComment("[pTError/pT]*max(1,normChi2) <= ptErrorCut");
  desc.add<std::string>("quality", "highPurity")->setComment(" quality cut as defined in reco::TrackBase");
  desc.add<bool>("useVtx", true)->setComment("compatibility with a vertex");
  desc.add<uint32_t>("nVertices", 0)->setComment(" how many vertices to look at before dropping the track");
  desc.add<edm::InputTag>("vertexTag", edm::InputTag("offlinePrimaryVertices"));
  desc.add<edm::InputTag>("timesTag", edm::InputTag(""));
  desc.add<edm::InputTag>("timeResosTag", edm::InputTag(""));
  desc.add<bool>("vtxFallback", true)->setComment("falback to beam spot if there are no vertices");
  desc.add<double>("zetaVtx", 1.0);
  desc.add<double>("rhoVtx", 0.2)->setComment("tags used by b-tagging folks");
  desc.add<double>("nSigmaDtVertex", 0);
}