TrajSeedMatcher.cc

Go to the documentation of this file.

#include "RecoEgamma/EgammaElectronAlgos/interface/TrajSeedMatcher.h"

#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"

#include "DataFormats/TrajectorySeed/interface/TrajectorySeed.h"
#include "DataFormats/GeometryCommonDetAlgo/interface/PerpendicularBoundPlaneBuilder.h"
#include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"

#include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
#include "TrackingTools/RecoGeometry/interface/RecoGeometryRecord.h"
#include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"

#include "RecoTracker/Record/interface/NavigationSchoolRecord.h"

#include "RecoEgamma/EgammaElectronAlgos/interface/FTSFromVertexToPointFactory.h"
#include "RecoEgamma/EgammaElectronAlgos/interface/ElectronUtilities.h"

constexpr float TrajSeedMatcher::kElectronMass_;

TrajSeedMatcher::TrajSeedMatcher(const edm::ParameterSet& pset)
    : cacheIDMagField_(0),
      paramMagFieldLabel_(pset.getParameter<std::string>("paramMagField")),
      navSchoolLabel_(pset.getParameter<std::string>("navSchool")),
      detLayerGeomLabel_(pset.getParameter<std::string>("detLayerGeom")),
      useRecoVertex_(pset.getParameter<bool>("useRecoVertex")),
      enableHitSkipping_(pset.getParameter<bool>("enableHitSkipping")),
      requireExactMatchCount_(pset.getParameter<bool>("requireExactMatchCount")),
      useParamMagFieldIfDefined_(pset.getParameter<bool>("useParamMagFieldIfDefined")),
      minNrHits_(pset.getParameter<std::vector<unsigned int> >("minNrHits")),
      minNrHitsValidLayerBins_(pset.getParameter<std::vector<int> >("minNrHitsValidLayerBins")) {
  const auto cutsPSets = pset.getParameter<std::vector<edm::ParameterSet> >("matchingCuts");
  for (const auto& cutPSet : cutsPSets) {
    int version = cutPSet.getParameter<int>("version");
    switch (version) {
      case 1:
        matchingCuts_.emplace_back(std::make_unique<MatchingCutsV1>(cutPSet));
        break;
      case 2:
        matchingCuts_.emplace_back(std::make_unique<MatchingCutsV2>(cutPSet));
        break;
      default:
        throw cms::Exception("InvalidConfig") << " Error TrajSeedMatcher::TrajSeedMatcher pixel match cuts version "
                                              << version << " not recognised" << std::endl;
    }
  }

  if (minNrHitsValidLayerBins_.size() + 1 != minNrHits_.size()) {
    throw cms::Exception("InvalidConfig")
        << " TrajSeedMatcher::TrajSeedMatcher minNrHitsValidLayerBins should be 1 less than minNrHits when its "
        << minNrHitsValidLayerBins_.size() << " vs " << minNrHits_.size();
  }
}

edm::ParameterSetDescription TrajSeedMatcher::makePSetDescription() {
  edm::ParameterSetDescription desc;
  desc.add<bool>("useRecoVertex", false);
  desc.add<bool>("enableHitSkipping", false);
  desc.add<bool>("requireExactMatchCount", true);
  desc.add<bool>("useParamMagFieldIfDefined", true);
  desc.add<std::string>("paramMagField", "ParabolicMf");
  desc.add<std::string>("navSchool", "SimpleNavigationSchool");
  desc.add<std::string>("detLayerGeom", "hltESPGlobalDetLayerGeometry");
  desc.add<std::vector<int> >("minNrHitsValidLayerBins", {4});
  desc.add<std::vector<unsigned int> >("minNrHits", {2, 3});

  edm::ParameterSetDescription cutsDesc;
  auto cutDescCases = 1 >> (edm::ParameterDescription<double>("dPhiMax", 0.04, true) and
                            edm::ParameterDescription<double>("dRZMax", 0.09, true) and
                            edm::ParameterDescription<double>("dRZMaxLowEtThres", 20., true) and
                            edm::ParameterDescription<std::vector<double> >(
                                "dRZMaxLowEtEtaBins", std::vector<double>{1., 1.5}, true) and
                            edm::ParameterDescription<std::vector<double> >(
                                "dRZMaxLowEt", std::vector<double>{0.09, 0.15, 0.09}, true)) or
                      2 >> (edm::ParameterDescription<std::vector<double> >("dPhiMaxHighEt", {0.003}, true) and
                            edm::ParameterDescription<std::vector<double> >("dPhiMaxHighEtThres", {0.0}, true) and
                            edm::ParameterDescription<std::vector<double> >("dPhiMaxLowEtGrad", {0.0}, true) and
                            edm::ParameterDescription<std::vector<double> >("dRZMaxHighEt", {0.005}, true) and
                            edm::ParameterDescription<std::vector<double> >("dRZMaxHighEtThres", {30}, true) and
                            edm::ParameterDescription<std::vector<double> >("dRZMaxLowEtGrad", {-0.002}, true) and
                            edm::ParameterDescription<std::vector<double> >("etaBins", {}, true));
  cutsDesc.ifValue(edm::ParameterDescription<int>("version", 1, true), std::move(cutDescCases));

  edm::ParameterSet defaults;
  defaults.addParameter<double>("dPhiMax", 0.04);
  defaults.addParameter<double>("dRZMax", 0.09);
  defaults.addParameter<double>("dRZMaxLowEtThres", 0.09);
  defaults.addParameter<std::vector<double> >("dRZMaxLowEtEtaBins", std::vector<double>{1., 1.5});
  defaults.addParameter<std::vector<double> >("dRZMaxLowEt", std::vector<double>{0.09, 0.09, 0.09});
  defaults.addParameter<int>("version", 1);
  desc.addVPSet("matchingCuts", cutsDesc, std::vector<edm::ParameterSet>{defaults, defaults, defaults});
  return desc;
}

void TrajSeedMatcher::doEventSetup(const edm::EventSetup& iSetup) {
  if (cacheIDMagField_ != iSetup.get<IdealMagneticFieldRecord>().cacheIdentifier()) {
    iSetup.get<IdealMagneticFieldRecord>().get(magField_);
    if (useParamMagFieldIfDefined_)
      iSetup.get<IdealMagneticFieldRecord>().get(paramMagFieldLabel_, magFieldParam_);
    cacheIDMagField_ = iSetup.get<IdealMagneticFieldRecord>().cacheIdentifier();
    forwardPropagator_ = std::make_unique<PropagatorWithMaterial>(alongMomentum, kElectronMass_, &*(magField_));
    backwardPropagator_ = std::make_unique<PropagatorWithMaterial>(oppositeToMomentum, kElectronMass_, &*(magField_));
  }
  iSetup.get<NavigationSchoolRecord>().get(navSchoolLabel_, navSchool_);
  iSetup.get<RecoGeometryRecord>().get(detLayerGeomLabel_, detLayerGeom_);
}

std::vector<TrajSeedMatcher::SeedWithInfo> TrajSeedMatcher::compatibleSeeds(const TrajectorySeedCollection& seeds,
                                                                            const GlobalPoint& candPos,
                                                                            const GlobalPoint& vprim,
                                                                            const float energy) {
  if (!forwardPropagator_ || !backwardPropagator_ || !magField_.isValid()) {
    throw cms::Exception("LogicError") << __FUNCTION__
                                       << " can not make pixel seeds as event setup has not properly been called";
  }

  clearCache();

  std::vector<SeedWithInfo> matchedSeeds;

  //these are super expensive functions
  TrajectoryStateOnSurface scTrajStateOnSurfNeg = makeTrajStateOnSurface(candPos, vprim, energy, -1);
  TrajectoryStateOnSurface scTrajStateOnSurfPos = makeTrajStateOnSurface(candPos, vprim, energy, 1);

  for (const auto& seed : seeds) {
    std::vector<SCHitMatch> matchedHitsNeg = processSeed(seed, candPos, vprim, energy, scTrajStateOnSurfNeg);
    std::vector<SCHitMatch> matchedHitsPos = processSeed(seed, candPos, vprim, energy, scTrajStateOnSurfPos);

    int nrValidLayersPos = 0;
    int nrValidLayersNeg = 0;
    if (matchedHitsNeg.size() >= 2) {
      nrValidLayersNeg = getNrValidLayersAlongTraj(matchedHitsNeg[0], matchedHitsNeg[1], candPos, vprim, energy, -1);
    }
    if (matchedHitsPos.size() >= 2) {
      nrValidLayersPos = getNrValidLayersAlongTraj(matchedHitsPos[0], matchedHitsPos[1], candPos, vprim, energy, +1);
    }

    int nrValidLayers = std::max(nrValidLayersNeg, nrValidLayersPos);
    size_t nrHitsRequired = getNrHitsRequired(nrValidLayers);
    bool matchCountPasses;
    if (requireExactMatchCount_) {
      // If the input seed collection is not cross-cleaned, an exact match is necessary to
      // prevent redundant seeds.
      matchCountPasses = matchedHitsNeg.size() == nrHitsRequired || matchedHitsPos.size() == nrHitsRequired;
    } else {
      matchCountPasses = matchedHitsNeg.size() >= nrHitsRequired || matchedHitsPos.size() >= nrHitsRequired;
    }
    if (matchCountPasses) {
      matchedSeeds.push_back({seed, matchedHitsPos, matchedHitsNeg, nrValidLayers});
    }
  }
  return matchedSeeds;
}

std::vector<TrajSeedMatcher::SCHitMatch> TrajSeedMatcher::processSeed(const TrajectorySeed& seed,
                                                                      const GlobalPoint& candPos,
                                                                      const GlobalPoint& vprim,
                                                                      const float energy,
                                                                      const TrajectoryStateOnSurface& initialTrajState) {
  //next try passing these variables in once...
  const float candEta = candPos.eta();
  const float candEt = energy * std::sin(candPos.theta());
  const int charge = initialTrajState.charge();

  std::vector<SCHitMatch> matchedHits;
  FreeTrajectoryState firstHitFreeTraj;
  GlobalPoint prevHitPos;
  GlobalPoint vertex;
  for (size_t hitNr = 0; hitNr < matchingCuts_.size() && hitNr < seed.nHits(); hitNr++) {
    bool matched = false;
    if (matchedHits.empty()) {  // First hit is special
      SCHitMatch firstHit = matchFirstHit(seed, initialTrajState, vprim, *backwardPropagator_);
      firstHit.setExtra(candEt, candEta, candPos.phi(), charge, 1);
      if (passesMatchSel(firstHit, 0)) {
        matched = true;
        matchedHits.push_back(firstHit);
        //now we can figure out the z vertex
        double zVertex = useRecoVertex_ ? vprim.z() : getZVtxFromExtrapolation(vprim, firstHit.hitPos(), candPos);
        vertex = GlobalPoint(vprim.x(), vprim.y(), zVertex);
        firstHitFreeTraj =
            FTSFromVertexToPointFactory::get(getMagField(firstHit.hitPos()), firstHit.hitPos(), vertex, energy, charge);
        prevHitPos = firstHit.hitPos();
      }
    } else {  // All subsequent hits are handled the same
      SCHitMatch hit = match2ndToNthHit(seed, firstHitFreeTraj, hitNr, prevHitPos, vertex, *forwardPropagator_);
      hit.setExtra(candEt, candEta, candPos.phi(), charge, 1);
      if (passesMatchSel(hit, matchedHits.size())) {
        matched = true;
        matchedHits.push_back(hit);
        prevHitPos = hit.hitPos();
      }
    }
    if (!matched and !enableHitSkipping_)
      break;
  }
  return matchedHits;
}

// compute the z vertex from the candidate position and the found pixel hit
float TrajSeedMatcher::getZVtxFromExtrapolation(const GlobalPoint& primeVtxPos,
                                                const GlobalPoint& hitPos,
                                                const GlobalPoint& candPos) {
  auto sq = [](float x) { return x * x; };
  auto calRDiff = [sq](const GlobalPoint& p1, const GlobalPoint& p2) {
    return std::sqrt(sq(p2.x() - p1.x()) + sq(p2.y() - p1.y()));
  };
  const double r1Diff = calRDiff(primeVtxPos, hitPos);
  const double r2Diff = calRDiff(hitPos, candPos);
  return hitPos.z() - r1Diff * (candPos.z() - hitPos.z()) / r2Diff;
}

bool TrajSeedMatcher::passTrajPreSel(const GlobalPoint& hitPos, const GlobalPoint& candPos) const {
  float dt = hitPos.x() * candPos.x() + hitPos.y() * candPos.y();
  if (dt < 0)
    return false;
  if (dt < kPhiCut_ * (candPos.perp() * hitPos.perp()))
    return false;
  return true;
}

const TrajectoryStateOnSurface& TrajSeedMatcher::getTrajStateFromVtx(const TrackingRecHit& hit,
                                                                     const TrajectoryStateOnSurface& initialState,
                                                                     const PropagatorWithMaterial& propagator) {
  auto& trajStateFromVtxCache =
      initialState.charge() == 1 ? trajStateFromVtxPosChargeCache_ : trajStateFromVtxNegChargeCache_;

  auto key = hit.det()->gdetIndex();
  auto res = trajStateFromVtxCache.find(key);
  if (res != trajStateFromVtxCache.end())
    return res->second;
  else {  //doesnt exist, need to make it
    //FIXME: check for efficiency
    auto val = trajStateFromVtxCache.emplace(key, propagator.propagate(initialState, hit.det()->surface()));
    return val.first->second;
  }
}

const TrajectoryStateOnSurface& TrajSeedMatcher::getTrajStateFromPoint(const TrackingRecHit& hit,
                                                                       const FreeTrajectoryState& initialState,
                                                                       const GlobalPoint& point,
                                                                       const PropagatorWithMaterial& propagator) {
  auto& trajStateFromPointCache =
      initialState.charge() == 1 ? trajStateFromPointPosChargeCache_ : trajStateFromPointNegChargeCache_;

  auto key = std::make_pair(hit.det()->gdetIndex(), point);
  auto res = trajStateFromPointCache.find(key);
  if (res != trajStateFromPointCache.end())
    return res->second;
  else {  //doesnt exist, need to make it
    //FIXME: check for efficiency
    auto val = trajStateFromPointCache.emplace(key, propagator.propagate(initialState, hit.det()->surface()));
    return val.first->second;
  }
}

TrajectoryStateOnSurface TrajSeedMatcher::makeTrajStateOnSurface(const GlobalPoint& pos,
                                                                 const GlobalPoint& vtx,
                                                                 const float energy,
                                                                 const int charge) const {
  FreeTrajectoryState freeTS = FTSFromVertexToPointFactory::get(getMagField(pos), pos, vtx, energy, charge);
  PerpendicularBoundPlaneBuilder bpb;
  return TrajectoryStateOnSurface(freeTS, *bpb(freeTS.position(), freeTS.momentum()));
}

TrajSeedMatcher::SCHitMatch TrajSeedMatcher::matchFirstHit(const TrajectorySeed& seed,
                                                           const TrajectoryStateOnSurface& initialState,
                                                           const GlobalPoint& vtxPos,
                                                           const PropagatorWithMaterial& propagator) {
  const TrajectorySeed::range& hits = seed.recHits();
  auto hitIt = hits.first;

  if (hitIt->isValid()) {
    const TrajectoryStateOnSurface& trajStateFromVtx = getTrajStateFromVtx(*hitIt, initialState, propagator);
    if (trajStateFromVtx.isValid())
      return SCHitMatch(vtxPos, trajStateFromVtx, *hitIt);
  }
  return SCHitMatch();
}

TrajSeedMatcher::SCHitMatch TrajSeedMatcher::match2ndToNthHit(const TrajectorySeed& seed,
                                                              const FreeTrajectoryState& initialState,
                                                              const size_t hitNr,
                                                              const GlobalPoint& prevHitPos,
                                                              const GlobalPoint& vtxPos,
                                                              const PropagatorWithMaterial& propagator) {
  const TrajectorySeed::range& hits = seed.recHits();
  auto hitIt = hits.first + hitNr;

  if (hitIt->isValid()) {
    const TrajectoryStateOnSurface& trajState = getTrajStateFromPoint(*hitIt, initialState, prevHitPos, propagator);

    if (trajState.isValid()) {
      return SCHitMatch(vtxPos, trajState, *hitIt);
    }
  }
  return SCHitMatch();
}

void TrajSeedMatcher::clearCache() {
  trajStateFromVtxPosChargeCache_.clear();
  trajStateFromVtxNegChargeCache_.clear();
  trajStateFromPointPosChargeCache_.clear();
  trajStateFromPointNegChargeCache_.clear();
}

bool TrajSeedMatcher::passesMatchSel(const TrajSeedMatcher::SCHitMatch& hit, const size_t hitNr) const {
  if (hitNr < matchingCuts_.size()) {
    return (*matchingCuts_[hitNr])(hit);
  } else {
    throw cms::Exception("LogicError") << " Error, attempting to apply selection to hit " << hitNr
                                       << " but only cuts for " << matchingCuts_.size() << " defined";
  }
}

int TrajSeedMatcher::getNrValidLayersAlongTraj(const SCHitMatch& hit1,
                                               const SCHitMatch& hit2,
                                               const GlobalPoint& candPos,
                                               const GlobalPoint& vprim,
                                               const float energy,
                                               const int charge) {
  double zVertex = useRecoVertex_ ? vprim.z() : getZVtxFromExtrapolation(vprim, hit1.hitPos(), candPos);
  GlobalPoint vertex(vprim.x(), vprim.y(), zVertex);

  FreeTrajectoryState firstHitFreeTraj =
      FTSFromVertexToPointFactory::get(getMagField(hit1.hitPos()), hit1.hitPos(), vertex, energy, charge);
  const TrajectoryStateOnSurface& secondHitTraj =
      getTrajStateFromPoint(*hit2.hit(), firstHitFreeTraj, hit1.hitPos(), *forwardPropagator_);
  return getNrValidLayersAlongTraj(hit2.hit()->geographicalId(), secondHitTraj);
}

int TrajSeedMatcher::getNrValidLayersAlongTraj(const DetId& hitId, const TrajectoryStateOnSurface& hitTrajState) const {
  const DetLayer* detLayer = detLayerGeom_->idToLayer(hitId);
  if (detLayer == nullptr)
    return 0;

  const FreeTrajectoryState& hitFreeState = *hitTrajState.freeState();
  const std::vector<const DetLayer*> inLayers =
      navSchool_->compatibleLayers(*detLayer, hitFreeState, oppositeToMomentum);
  const std::vector<const DetLayer*> outLayers = navSchool_->compatibleLayers(*detLayer, hitFreeState, alongMomentum);

  int nrValidLayers = 1;  //because our current hit is also valid and wont be included in the count otherwise
  int nrPixInLayers = 0;
  int nrPixOutLayers = 0;
  for (auto layer : inLayers) {
    if (GeomDetEnumerators::isTrackerPixel(layer->subDetector())) {
      nrPixInLayers++;
      if (layerHasValidHits(*layer, hitTrajState, *backwardPropagator_))
        nrValidLayers++;
    }
  }
  for (auto layer : outLayers) {
    if (GeomDetEnumerators::isTrackerPixel(layer->subDetector())) {
      nrPixOutLayers++;
      if (layerHasValidHits(*layer, hitTrajState, *forwardPropagator_))
        nrValidLayers++;
    }
  }
  return nrValidLayers;
}

bool TrajSeedMatcher::layerHasValidHits(const DetLayer& layer,
                                        const TrajectoryStateOnSurface& hitSurState,
                                        const Propagator& propToLayerFromState) const {
  //FIXME: do not hardcode with werid magic numbers stolen from ancient tracking code
  //its taken from https://cmssdt.cern.ch/dxr/CMSSW/source/RecoTracker/TrackProducer/interface/TrackProducerBase.icc#165
  //which inspires this code
  Chi2MeasurementEstimator estimator(30., -3.0, 0.5, 2.0, 0.5, 1.e12);  // same as defauts....

  const std::vector<GeometricSearchDet::DetWithState>& detWithState =
      layer.compatibleDets(hitSurState, propToLayerFromState, estimator);
  if (detWithState.empty())
    return false;
  else {
    DetId id = detWithState.front().first->geographicalId();
    MeasurementDetWithData measDet = measTkEvt_->idToDet(id);
    if (measDet.isActive())
      return true;
    else
      return false;
  }
}

size_t TrajSeedMatcher::getNrHitsRequired(const int nrValidLayers) const {
  for (size_t binNr = 0; binNr < minNrHitsValidLayerBins_.size(); binNr++) {
    if (nrValidLayers < minNrHitsValidLayerBins_[binNr])
      return minNrHits_[binNr];
  }
  return minNrHits_.back();
}

TrajSeedMatcher::SCHitMatch::SCHitMatch(const GlobalPoint& vtxPos,
                                        const TrajectoryStateOnSurface& trajState,
                                        const TrackingRecHit& hit)
    : detId_(hit.geographicalId()),
      hitPos_(hit.globalPosition()),
      hit_(&hit),
      et_(0),
      eta_(0),
      phi_(0),
      charge_(0),
      nrClus_(0) {
  EleRelPointPair pointPair(hitPos_, trajState.globalParameters().position(), vtxPos);
  dRZ_ = detId_.subdetId() == PixelSubdetector::PixelBarrel ? pointPair.dZ() : pointPair.dPerp();
  dPhi_ = pointPair.dPhi();
}

TrajSeedMatcher::SeedWithInfo::SeedWithInfo(const TrajectorySeed& seed,
                                            const std::vector<SCHitMatch>& posCharge,
                                            const std::vector<SCHitMatch>& negCharge,
                                            int nrValidLayers)
    : seed_(seed), nrValidLayers_(nrValidLayers) {
  size_t nrHitsMax = std::max(posCharge.size(), negCharge.size());
  for (size_t hitNr = 0; hitNr < nrHitsMax; hitNr++) {
    DetId detIdPos = hitNr < posCharge.size() ? posCharge[hitNr].detId() : DetId(0);
    float dRZPos = hitNr < posCharge.size() ? posCharge[hitNr].dRZ() : std::numeric_limits<float>::max();
    float dPhiPos = hitNr < posCharge.size() ? posCharge[hitNr].dPhi() : std::numeric_limits<float>::max();

    DetId detIdNeg = hitNr < negCharge.size() ? negCharge[hitNr].detId() : DetId(0);
    float dRZNeg = hitNr < negCharge.size() ? negCharge[hitNr].dRZ() : std::numeric_limits<float>::max();
    float dPhiNeg = hitNr < negCharge.size() ? negCharge[hitNr].dPhi() : std::numeric_limits<float>::max();

    if (detIdPos != detIdNeg && (detIdPos.rawId() != 0 && detIdNeg.rawId() != 0)) {
      cms::Exception("LogicError")
          << " error in " << __FILE__ << ", " << __LINE__
          << " hits to be combined have different detIDs, this should not be possible and nothing good will come of it";
    }
    DetId detId = detIdPos.rawId() != 0 ? detIdPos : detIdNeg;
    matchInfo_.push_back(MatchInfo(detId, dRZPos, dRZNeg, dPhiPos, dPhiNeg));
  }
}

TrajSeedMatcher::MatchingCutsV1::MatchingCutsV1(const edm::ParameterSet& pset)
    : dPhiMax_(pset.getParameter<double>("dPhiMax")),
      dRZMax_(pset.getParameter<double>("dRZMax")),
      dRZMaxLowEtThres_(pset.getParameter<double>("dRZMaxLowEtThres")),
      dRZMaxLowEtEtaBins_(pset.getParameter<std::vector<double> >("dRZMaxLowEtEtaBins")),
      dRZMaxLowEt_(pset.getParameter<std::vector<double> >("dRZMaxLowEt")) {
  if (dRZMaxLowEtEtaBins_.size() + 1 != dRZMaxLowEt_.size()) {
    throw cms::Exception("InvalidConfig") << " dRZMaxLowEtEtaBins should be 1 less than dRZMaxLowEt when its "
                                          << dRZMaxLowEtEtaBins_.size() << " vs " << dRZMaxLowEt_.size();
  }
}

bool TrajSeedMatcher::MatchingCutsV1::operator()(const TrajSeedMatcher::SCHitMatch& scHitMatch) const {
  if (dPhiMax_ >= 0 && std::abs(scHitMatch.dPhi()) > dPhiMax_)
    return false;

  const float dRZMax = getDRZCutValue(scHitMatch.et(), scHitMatch.eta());
  if (dRZMax_ >= 0 && std::abs(scHitMatch.dRZ()) > dRZMax)
    return false;

  return true;
}

float TrajSeedMatcher::MatchingCutsV1::getDRZCutValue(const float scEt, const float scEta) const {
  if (scEt >= dRZMaxLowEtThres_)
    return dRZMax_;
  else {
    const float absEta = std::abs(scEta);
    for (size_t etaNr = 0; etaNr < dRZMaxLowEtEtaBins_.size(); etaNr++) {
      if (absEta < dRZMaxLowEtEtaBins_[etaNr])
        return dRZMaxLowEt_[etaNr];
    }
    return dRZMaxLowEt_.back();
  }
}

TrajSeedMatcher::MatchingCutsV2::MatchingCutsV2(const edm::ParameterSet& pset)
    : dPhiHighEt_(pset.getParameter<std::vector<double> >("dPhiMaxHighEt")),
      dPhiHighEtThres_(pset.getParameter<std::vector<double> >("dPhiMaxHighEtThres")),
      dPhiLowEtGrad_(pset.getParameter<std::vector<double> >("dPhiMaxLowEtGrad")),
      dRZHighEt_(pset.getParameter<std::vector<double> >("dRZMaxHighEt")),
      dRZHighEtThres_(pset.getParameter<std::vector<double> >("dRZMaxHighEtThres")),
      dRZLowEtGrad_(pset.getParameter<std::vector<double> >("dRZMaxLowEtGrad")),
      etaBins_(pset.getParameter<std::vector<double> >("etaBins")) {
  auto binSizeCheck = [](size_t sizeEtaBins, const std::vector<double>& vec, const std::string& name) {
    if (vec.size() != sizeEtaBins + 1) {
      throw cms::Exception("InvalidConfig")
          << " when constructing TrajSeedMatcher::MatchingCutsV2 " << name << " has " << vec.size()
          << " bins, it should be equal to #bins of etaBins+1" << sizeEtaBins + 1;
    }
  };
  binSizeCheck(etaBins_.size(), dPhiHighEt_, "dPhiMaxHighEt");
  binSizeCheck(etaBins_.size(), dPhiHighEtThres_, "dPhiMaxHighEtThres");
  binSizeCheck(etaBins_.size(), dPhiLowEtGrad_, "dPhiMaxLowEtGrad");
  binSizeCheck(etaBins_.size(), dRZHighEt_, "dRZMaxHighEt");
  binSizeCheck(etaBins_.size(), dRZHighEtThres_, "dRZMaxHighEtThres");
  binSizeCheck(etaBins_.size(), dRZLowEtGrad_, "dRZMaxLowEtGrad");
}

bool TrajSeedMatcher::MatchingCutsV2::operator()(const TrajSeedMatcher::SCHitMatch& scHitMatch) const {
  size_t binNr = getBinNr(scHitMatch.eta());
  float dPhiMax = getCutValue(scHitMatch.et(), dPhiHighEt_[binNr], dPhiHighEtThres_[binNr], dPhiLowEtGrad_[binNr]);
  if (dPhiMax >= 0 && std::abs(scHitMatch.dPhi()) > dPhiMax)
    return false;
  float dRZMax = getCutValue(scHitMatch.et(), dRZHighEt_[binNr], dRZHighEtThres_[binNr], dRZLowEtGrad_[binNr]);
  if (dRZMax >= 0 && std::abs(scHitMatch.dRZ()) > dRZMax)
    return false;

  return true;
}

//eta bins is exactly 1 smaller than the vectors which will be accessed by this bin nr
size_t TrajSeedMatcher::MatchingCutsV2::getBinNr(float eta) const {
  const float absEta = std::abs(eta);
  for (size_t etaNr = 0; etaNr < etaBins_.size(); etaNr++) {
    if (absEta < etaBins_[etaNr])
      return etaNr;
  }
  return etaBins_.size();
}