Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:34:03

0001 #include "RecoEgamma/EgammaElectronAlgos/interface/TrajSeedMatcher.h"
0002 
0003 #include "FWCore/Framework/interface/EventSetup.h"
0004 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0005 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0006 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0007 
0008 #include "DataFormats/TrajectorySeed/interface/TrajectorySeed.h"
0009 #include "DataFormats/GeometryCommonDetAlgo/interface/PerpendicularBoundPlaneBuilder.h"
0010 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0011 
0012 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0013 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"
0014 
0015 #include "TrackingTools/TrajectoryState/interface/ftsFromVertexToPoint.h"
0016 #include "RecoEgamma/EgammaElectronAlgos/interface/ElectronUtilities.h"
0017 
0018 #include "FWCore/Framework/interface/ConsumesCollector.h"
0019 
0020 constexpr float TrajSeedMatcher::kElectronMass_;
0021 
0022 namespace {
0023   auto makeMatchingCuts(std::vector<edm::ParameterSet> const& cutsPSets) {
0024     std::vector<std::unique_ptr<TrajSeedMatcher::MatchingCuts> > matchingCuts;
0025 
0026     for (const auto& cutPSet : cutsPSets) {
0027       int version = cutPSet.getParameter<int>("version");
0028       switch (version) {
0029         case 1:
0030           matchingCuts.emplace_back(std::make_unique<TrajSeedMatcher::MatchingCutsV1>(cutPSet));
0031           break;
0032         case 2:
0033           matchingCuts.emplace_back(std::make_unique<TrajSeedMatcher::MatchingCutsV2>(cutPSet));
0034           break;
0035         default:
0036           throw cms::Exception("InvalidConfig") << " Error TrajSeedMatcher::TrajSeedMatcher pixel match cuts version "
0037                                                 << version << " not recognised" << std::endl;
0038       }
0039     }
0040 
0041     return matchingCuts;
0042   }
0043 
0044   TrajSeedMatcher::SCHitMatch makeSCHitMatch(const GlobalPoint& vtxPos,
0045                                              const TrajectoryStateOnSurface& trajState,
0046                                              const TrackingRecHit& hit,
0047                                              float et,
0048                                              float eta,
0049                                              float phi,
0050                                              int charge,
0051                                              int nrClus) {
0052     EleRelPointPair pointPair(hit.globalPosition(), trajState.globalParameters().position(), vtxPos);
0053     float dRZ = hit.geographicalId().subdetId() == PixelSubdetector::PixelBarrel ? pointPair.dZ() : pointPair.dPerp();
0054     return {hit.geographicalId(), hit.globalPosition(), dRZ, pointPair.dPhi(), hit, et, eta, phi, charge, nrClus};
0055   }
0056 
0057   const std::vector<TrajSeedMatcher::MatchInfo> makeMatchInfoVector(
0058       std::vector<TrajSeedMatcher::SCHitMatch> const& posCharge,
0059       std::vector<TrajSeedMatcher::SCHitMatch> const& negCharge) {
0060     std::vector<TrajSeedMatcher::MatchInfo> matchInfos;
0061     size_t nrHitsMax = std::max(posCharge.size(), negCharge.size());
0062     for (size_t hitNr = 0; hitNr < nrHitsMax; hitNr++) {
0063       DetId detIdPos = hitNr < posCharge.size() ? posCharge[hitNr].detId : DetId(0);
0064       float dRZPos = hitNr < posCharge.size() ? posCharge[hitNr].dRZ : std::numeric_limits<float>::max();
0065       float dPhiPos = hitNr < posCharge.size() ? posCharge[hitNr].dPhi : std::numeric_limits<float>::max();
0066 
0067       DetId detIdNeg = hitNr < negCharge.size() ? negCharge[hitNr].detId : DetId(0);
0068       float dRZNeg = hitNr < negCharge.size() ? negCharge[hitNr].dRZ : std::numeric_limits<float>::max();
0069       float dPhiNeg = hitNr < negCharge.size() ? negCharge[hitNr].dPhi : std::numeric_limits<float>::max();
0070 
0071       if (detIdPos != detIdNeg && (detIdPos.rawId() != 0 && detIdNeg.rawId() != 0)) {
0072         cms::Exception("LogicError") << " error in " << __FILE__ << ", " << __LINE__
0073                                      << " hits to be combined have different detIDs, this should not be possible and "
0074                                         "nothing good will come of it";
0075       }
0076       DetId detId = detIdPos.rawId() != 0 ? detIdPos : detIdNeg;
0077       matchInfos.push_back({detId, dRZPos, dRZNeg, dPhiPos, dPhiNeg});
0078     }
0079     return matchInfos;
0080   }
0081 };  // namespace
0082 
0083 TrajSeedMatcher::Configuration::Configuration(const edm::ParameterSet& pset, edm::ConsumesCollector&& cc)
0084     : magFieldToken{cc.esConsumes()},
0085       paramMagFieldToken{cc.esConsumes(pset.getParameter<edm::ESInputTag>("paramMagField"))},
0086       navSchoolToken{cc.esConsumes(pset.getParameter<edm::ESInputTag>("navSchool"))},
0087       detLayerGeomToken{cc.esConsumes(pset.getParameter<edm::ESInputTag>("detLayerGeom"))},
0088       useRecoVertex{pset.getParameter<bool>("useRecoVertex")},
0089       enableHitSkipping{pset.getParameter<bool>("enableHitSkipping")},
0090       requireExactMatchCount{pset.getParameter<bool>("requireExactMatchCount")},
0091       useParamMagFieldIfDefined{pset.getParameter<bool>("useParamMagFieldIfDefined")},
0092       minNrHits{pset.getParameter<std::vector<unsigned int> >("minNrHits")},
0093       minNrHitsValidLayerBins{pset.getParameter<std::vector<int> >("minNrHitsValidLayerBins")},
0094       matchingCuts{makeMatchingCuts(pset.getParameter<std::vector<edm::ParameterSet> >("matchingCuts"))} {
0095   if (minNrHitsValidLayerBins.size() + 1 != minNrHits.size()) {
0096     throw cms::Exception("InvalidConfig")
0097         << " TrajSeedMatcher::TrajSeedMatcher minNrHitsValidLayerBins should be 1 less than minNrHits when its "
0098         << minNrHitsValidLayerBins.size() << " vs " << minNrHits.size();
0099   }
0100 }
0101 
0102 TrajSeedMatcher::TrajSeedMatcher(TrajectorySeedCollection const& seeds,
0103                                  math::XYZPoint const& vprim,
0104                                  TrajSeedMatcher::Configuration const& cfg,
0105                                  edm::EventSetup const& iSetup,
0106                                  MeasurementTrackerEvent const& measTkEvt)
0107     : seeds_{seeds},
0108       vprim_(vprim.x(), vprim.y(), vprim.z()),
0109       cfg_{cfg},
0110       magField_{iSetup.getData(cfg_.magFieldToken)},
0111       magFieldParam_{iSetup.getData(cfg_.paramMagFieldToken)},
0112       measTkEvt_{measTkEvt},
0113       navSchool_{iSetup.getData(cfg_.navSchoolToken)},
0114       detLayerGeom_{iSetup.getData(cfg_.detLayerGeomToken)},
0115       forwardPropagator_(alongMomentum, kElectronMass_, &magField_),
0116       backwardPropagator_(oppositeToMomentum, kElectronMass_, &magField_) {}
0117 
0118 edm::ParameterSetDescription TrajSeedMatcher::makePSetDescription() {
0119   edm::ParameterSetDescription desc;
0120   desc.add<bool>("useRecoVertex", false);
0121   desc.add<bool>("enableHitSkipping", false);
0122   desc.add<bool>("requireExactMatchCount", true);
0123   desc.add<bool>("useParamMagFieldIfDefined", true);
0124   desc.add<edm::ESInputTag>("paramMagField", edm::ESInputTag{"", "ParabolicMf"});
0125   desc.add<edm::ESInputTag>("navSchool", edm::ESInputTag{"", "SimpleNavigationSchool"});
0126   desc.add<edm::ESInputTag>("detLayerGeom", edm::ESInputTag{"", "hltESPGlobalDetLayerGeometry"});
0127   desc.add<std::vector<int> >("minNrHitsValidLayerBins", {4});
0128   desc.add<std::vector<unsigned int> >("minNrHits", {2, 3});
0129 
0130   edm::ParameterSetDescription cutsDesc;
0131   auto cutDescCases = 1 >> (edm::ParameterDescription<double>("dPhiMax", 0.04, true) and
0132                             edm::ParameterDescription<double>("dRZMax", 0.09, true) and
0133                             edm::ParameterDescription<double>("dRZMaxLowEtThres", 20., true) and
0134                             edm::ParameterDescription<std::vector<double> >("dRZMaxLowEtEtaBins", {1., 1.5}, true) and
0135                             edm::ParameterDescription<std::vector<double> >("dRZMaxLowEt", {0.09, 0.15, 0.09}, true)) or
0136                       2 >> (edm::ParameterDescription<std::vector<double> >("dPhiMaxHighEt", {0.003}, true) and
0137                             edm::ParameterDescription<std::vector<double> >("dPhiMaxHighEtThres", {0.0}, true) and
0138                             edm::ParameterDescription<std::vector<double> >("dPhiMaxLowEtGrad", {0.0}, true) and
0139                             edm::ParameterDescription<std::vector<double> >("dRZMaxHighEt", {0.005}, true) and
0140                             edm::ParameterDescription<std::vector<double> >("dRZMaxHighEtThres", {30}, true) and
0141                             edm::ParameterDescription<std::vector<double> >("dRZMaxLowEtGrad", {-0.002}, true) and
0142                             edm::ParameterDescription<std::vector<double> >("etaBins", {}, true));
0143   cutsDesc.ifValue(edm::ParameterDescription<int>("version", 1, true), std::move(cutDescCases));
0144 
0145   edm::ParameterSet defaults;
0146   defaults.addParameter<double>("dPhiMax", 0.04);
0147   defaults.addParameter<double>("dRZMax", 0.09);
0148   defaults.addParameter<double>("dRZMaxLowEtThres", 0.09);
0149   defaults.addParameter<std::vector<double> >("dRZMaxLowEtEtaBins", std::vector<double>{1., 1.5});
0150   defaults.addParameter<std::vector<double> >("dRZMaxLowEt", std::vector<double>{0.09, 0.09, 0.09});
0151   defaults.addParameter<int>("version", 1);
0152   desc.addVPSet("matchingCuts", cutsDesc, std::vector<edm::ParameterSet>{defaults, defaults, defaults});
0153   return desc;
0154 }
0155 
0156 std::vector<TrajSeedMatcher::SeedWithInfo> TrajSeedMatcher::operator()(const GlobalPoint& candPos, const float energy) {
0157   clearCache();
0158 
0159   std::vector<SeedWithInfo> matchedSeeds;
0160 
0161   //these are super expensive functions
0162   TrajectoryStateOnSurface scTrajStateOnSurfNeg = makeTrajStateOnSurface(candPos, energy, -1);
0163   TrajectoryStateOnSurface scTrajStateOnSurfPos = makeTrajStateOnSurface(candPos, energy, 1);
0164 
0165   for (const auto& seed : seeds_) {
0166     std::vector<SCHitMatch> matchedHitsNeg = processSeed(seed, candPos, energy, scTrajStateOnSurfNeg);
0167     std::vector<SCHitMatch> matchedHitsPos = processSeed(seed, candPos, energy, scTrajStateOnSurfPos);
0168 
0169     int nrValidLayersPos = 0;
0170     int nrValidLayersNeg = 0;
0171     if (matchedHitsNeg.size() >= 2) {
0172       nrValidLayersNeg = getNrValidLayersAlongTraj(matchedHitsNeg[0], matchedHitsNeg[1], candPos, energy, -1);
0173     }
0174     if (matchedHitsPos.size() >= 2) {
0175       nrValidLayersPos = getNrValidLayersAlongTraj(matchedHitsPos[0], matchedHitsPos[1], candPos, energy, +1);
0176     }
0177 
0178     int nrValidLayers = std::max(nrValidLayersNeg, nrValidLayersPos);
0179     size_t nrHitsRequired = getNrHitsRequired(nrValidLayers);
0180     bool matchCountPasses;
0181     if (cfg_.requireExactMatchCount) {
0182       // If the input seed collection is not cross-cleaned, an exact match is necessary to
0183       // prevent redundant seeds.
0184       matchCountPasses = matchedHitsNeg.size() == nrHitsRequired || matchedHitsPos.size() == nrHitsRequired;
0185     } else {
0186       matchCountPasses = matchedHitsNeg.size() >= nrHitsRequired || matchedHitsPos.size() >= nrHitsRequired;
0187     }
0188     if (matchCountPasses) {
0189       matchedSeeds.push_back({seed, makeMatchInfoVector(matchedHitsPos, matchedHitsNeg), nrValidLayers});
0190     }
0191   }
0192   return matchedSeeds;
0193 }
0194 
0195 std::vector<TrajSeedMatcher::SCHitMatch> TrajSeedMatcher::processSeed(const TrajectorySeed& seed,
0196                                                                       const GlobalPoint& candPos,
0197                                                                       const float energy,
0198                                                                       const TrajectoryStateOnSurface& initialTrajState) {
0199   //next try passing these variables in once...
0200   const float candEta = candPos.eta();
0201   const float candEt = energy * std::sin(candPos.theta());
0202   const int charge = initialTrajState.charge();
0203 
0204   std::vector<SCHitMatch> matches;
0205   FreeTrajectoryState firstMatchFreeTraj;
0206   GlobalPoint prevHitPos;
0207   GlobalPoint vertex;
0208   const auto nCuts = cfg_.matchingCuts.size();
0209   for (size_t iHit = 0;
0210        matches.size() < nCuts && iHit < seed.nHits() && (cfg_.enableHitSkipping || iHit == matches.size());
0211        iHit++) {
0212     auto const& recHit = *(seed.recHits().begin() + iHit);
0213 
0214     if (!recHit.isValid()) {
0215       continue;
0216     }
0217 
0218     const bool doFirstMatch = matches.empty();
0219 
0220     auto const& trajState = doFirstMatch
0221                                 ? getTrajStateFromVtx(recHit, initialTrajState, backwardPropagator_)
0222                                 : getTrajStateFromPoint(recHit, firstMatchFreeTraj, prevHitPos, forwardPropagator_);
0223     if (!trajState.isValid()) {
0224       continue;
0225     }
0226 
0227     auto const& vtxForMatchObject = doFirstMatch ? vprim_ : vertex;
0228     auto match = makeSCHitMatch(vtxForMatchObject, trajState, recHit, candEt, candEta, candPos.phi(), charge, 1);
0229 
0230     if ((*cfg_.matchingCuts[matches.size()])(match)) {
0231       matches.push_back(match);
0232       if (doFirstMatch) {
0233         //now we can figure out the z vertex
0234         double zVertex = cfg_.useRecoVertex ? vprim_.z() : getZVtxFromExtrapolation(vprim_, match.hitPos, candPos);
0235         vertex = GlobalPoint(vprim_.x(), vprim_.y(), zVertex);
0236         firstMatchFreeTraj = ftsFromVertexToPoint(match.hitPos, vertex, energy, charge);
0237       }
0238       prevHitPos = match.hitPos;
0239     }
0240   }
0241   return matches;
0242 }
0243 
0244 // compute the z vertex from the candidate position and the found pixel hit
0245 float TrajSeedMatcher::getZVtxFromExtrapolation(const GlobalPoint& primeVtxPos,
0246                                                 const GlobalPoint& hitPos,
0247                                                 const GlobalPoint& candPos) {
0248   auto sq = [](float x) { return x * x; };
0249   auto calRDiff = [sq](const GlobalPoint& p1, const GlobalPoint& p2) {
0250     return std::sqrt(sq(p2.x() - p1.x()) + sq(p2.y() - p1.y()));
0251   };
0252   const double r1Diff = calRDiff(primeVtxPos, hitPos);
0253   const double r2Diff = calRDiff(hitPos, candPos);
0254   return hitPos.z() - r1Diff * (candPos.z() - hitPos.z()) / r2Diff;
0255 }
0256 
0257 const TrajectoryStateOnSurface& TrajSeedMatcher::getTrajStateFromVtx(const TrackingRecHit& hit,
0258                                                                      const TrajectoryStateOnSurface& initialState,
0259                                                                      const PropagatorWithMaterial& propagator) {
0260   auto& trajStateFromVtxCache =
0261       initialState.charge() == 1 ? trajStateFromVtxPosChargeCache_ : trajStateFromVtxNegChargeCache_;
0262 
0263   auto key = hit.det()->gdetIndex();
0264   auto res = trajStateFromVtxCache.find(key);
0265   if (res != trajStateFromVtxCache.end())
0266     return res->second;
0267   else {  //doesnt exist, need to make it
0268     //FIXME: check for efficiency
0269     auto val = trajStateFromVtxCache.emplace(key, propagator.propagate(initialState, hit.det()->surface()));
0270     return val.first->second;
0271   }
0272 }
0273 
0274 const TrajectoryStateOnSurface& TrajSeedMatcher::getTrajStateFromPoint(const TrackingRecHit& hit,
0275                                                                        const FreeTrajectoryState& initialState,
0276                                                                        const GlobalPoint& point,
0277                                                                        const PropagatorWithMaterial& propagator) {
0278   auto& trajStateFromPointCache =
0279       initialState.charge() == 1 ? trajStateFromPointPosChargeCache_ : trajStateFromPointNegChargeCache_;
0280 
0281   auto key = std::make_pair(hit.det()->gdetIndex(), point);
0282   auto res = trajStateFromPointCache.find(key);
0283   if (res != trajStateFromPointCache.end())
0284     return res->second;
0285   else {  //doesnt exist, need to make it
0286     //FIXME: check for efficiency
0287     auto val = trajStateFromPointCache.emplace(key, propagator.propagate(initialState, hit.det()->surface()));
0288     return val.first->second;
0289   }
0290 }
0291 
0292 TrajectoryStateOnSurface TrajSeedMatcher::makeTrajStateOnSurface(const GlobalPoint& pos,
0293                                                                  const float energy,
0294                                                                  const int charge) const {
0295   auto freeTS = ftsFromVertexToPoint(pos, vprim_, energy, charge);
0296   return TrajectoryStateOnSurface(freeTS, *PerpendicularBoundPlaneBuilder{}(freeTS.position(), freeTS.momentum()));
0297 }
0298 
0299 void TrajSeedMatcher::clearCache() {
0300   trajStateFromVtxPosChargeCache_.clear();
0301   trajStateFromVtxNegChargeCache_.clear();
0302   trajStateFromPointPosChargeCache_.clear();
0303   trajStateFromPointNegChargeCache_.clear();
0304 }
0305 
0306 int TrajSeedMatcher::getNrValidLayersAlongTraj(
0307     const SCHitMatch& hit1, const SCHitMatch& hit2, const GlobalPoint& candPos, const float energy, const int charge) {
0308   double zVertex = cfg_.useRecoVertex ? vprim_.z() : getZVtxFromExtrapolation(vprim_, hit1.hitPos, candPos);
0309   GlobalPoint vertex(vprim_.x(), vprim_.y(), zVertex);
0310 
0311   auto firstMatchFreeTraj = ftsFromVertexToPoint(hit1.hitPos, vertex, energy, charge);
0312   auto const& secondHitTraj = getTrajStateFromPoint(hit2.hit, firstMatchFreeTraj, hit1.hitPos, forwardPropagator_);
0313   return getNrValidLayersAlongTraj(hit2.hit.geographicalId(), secondHitTraj);
0314 }
0315 
0316 int TrajSeedMatcher::getNrValidLayersAlongTraj(const DetId& hitId, const TrajectoryStateOnSurface& hitTrajState) const {
0317   const DetLayer* detLayer = detLayerGeom_.idToLayer(hitId);
0318   if (detLayer == nullptr)
0319     return 0;
0320 
0321   const FreeTrajectoryState& hitFreeState = *hitTrajState.freeState();
0322   auto const inLayers = navSchool_.compatibleLayers(*detLayer, hitFreeState, oppositeToMomentum);
0323   const auto outLayers = navSchool_.compatibleLayers(*detLayer, hitFreeState, alongMomentum);
0324 
0325   int nrValidLayers = 1;  //because our current hit is also valid and wont be included in the count otherwise
0326   int nrPixInLayers = 0;
0327   int nrPixOutLayers = 0;
0328   for (auto layer : inLayers) {
0329     if (GeomDetEnumerators::isTrackerPixel(layer->subDetector())) {
0330       nrPixInLayers++;
0331       if (layerHasValidHits(*layer, hitTrajState, backwardPropagator_))
0332         nrValidLayers++;
0333     }
0334   }
0335   for (auto layer : outLayers) {
0336     if (GeomDetEnumerators::isTrackerPixel(layer->subDetector())) {
0337       nrPixOutLayers++;
0338       if (layerHasValidHits(*layer, hitTrajState, forwardPropagator_))
0339         nrValidLayers++;
0340     }
0341   }
0342   return nrValidLayers;
0343 }
0344 
0345 bool TrajSeedMatcher::layerHasValidHits(const DetLayer& layer,
0346                                         const TrajectoryStateOnSurface& hitSurState,
0347                                         const Propagator& propToLayerFromState) const {
0348   //FIXME: do not hardcode with werid magic numbers stolen from ancient tracking code
0349   //its taken from https://cmssdt.cern.ch/dxr/CMSSW/source/RecoTracker/TrackProducer/interface/TrackProducerBase.icc#165
0350   //which inspires this code
0351   Chi2MeasurementEstimator estimator(30., -3.0, 0.5, 2.0, 0.5, 1.e12);  // same as defauts....
0352 
0353   const std::vector<GeometricSearchDet::DetWithState>& detWithState =
0354       layer.compatibleDets(hitSurState, propToLayerFromState, estimator);
0355   if (detWithState.empty())
0356     return false;
0357   else {
0358     DetId id = detWithState.front().first->geographicalId();
0359     MeasurementDetWithData measDet = measTkEvt_.idToDet(id);
0360     if (measDet.isActive())
0361       return true;
0362     else
0363       return false;
0364   }
0365 }
0366 
0367 size_t TrajSeedMatcher::getNrHitsRequired(const int nrValidLayers) const {
0368   for (size_t binNr = 0; binNr < cfg_.minNrHitsValidLayerBins.size(); binNr++) {
0369     if (nrValidLayers < cfg_.minNrHitsValidLayerBins[binNr])
0370       return cfg_.minNrHits[binNr];
0371   }
0372   return cfg_.minNrHits.back();
0373 }
0374 
0375 TrajSeedMatcher::MatchingCutsV1::MatchingCutsV1(const edm::ParameterSet& pset)
0376     : dPhiMax_(pset.getParameter<double>("dPhiMax")),
0377       dRZMax_(pset.getParameter<double>("dRZMax")),
0378       dRZMaxLowEtThres_(pset.getParameter<double>("dRZMaxLowEtThres")),
0379       dRZMaxLowEtEtaBins_(pset.getParameter<std::vector<double> >("dRZMaxLowEtEtaBins")),
0380       dRZMaxLowEt_(pset.getParameter<std::vector<double> >("dRZMaxLowEt")) {
0381   if (dRZMaxLowEtEtaBins_.size() + 1 != dRZMaxLowEt_.size()) {
0382     throw cms::Exception("InvalidConfig") << " dRZMaxLowEtEtaBins should be 1 less than dRZMaxLowEt when its "
0383                                           << dRZMaxLowEtEtaBins_.size() << " vs " << dRZMaxLowEt_.size();
0384   }
0385 }
0386 
0387 bool TrajSeedMatcher::MatchingCutsV1::operator()(const TrajSeedMatcher::SCHitMatch& scHitMatch) const {
0388   if (dPhiMax_ >= 0 && std::abs(scHitMatch.dPhi) > dPhiMax_)
0389     return false;
0390 
0391   const float dRZMax = getDRZCutValue(scHitMatch.et, scHitMatch.eta);
0392   if (dRZMax_ >= 0 && std::abs(scHitMatch.dRZ) > dRZMax)
0393     return false;
0394 
0395   return true;
0396 }
0397 
0398 float TrajSeedMatcher::MatchingCutsV1::getDRZCutValue(const float scEt, const float scEta) const {
0399   if (scEt >= dRZMaxLowEtThres_)
0400     return dRZMax_;
0401   else {
0402     const float absEta = std::abs(scEta);
0403     for (size_t etaNr = 0; etaNr < dRZMaxLowEtEtaBins_.size(); etaNr++) {
0404       if (absEta < dRZMaxLowEtEtaBins_[etaNr])
0405         return dRZMaxLowEt_[etaNr];
0406     }
0407     return dRZMaxLowEt_.back();
0408   }
0409 }
0410 
0411 TrajSeedMatcher::MatchingCutsV2::MatchingCutsV2(const edm::ParameterSet& pset)
0412     : dPhiHighEt_(pset.getParameter<std::vector<double> >("dPhiMaxHighEt")),
0413       dPhiHighEtThres_(pset.getParameter<std::vector<double> >("dPhiMaxHighEtThres")),
0414       dPhiLowEtGrad_(pset.getParameter<std::vector<double> >("dPhiMaxLowEtGrad")),
0415       dRZHighEt_(pset.getParameter<std::vector<double> >("dRZMaxHighEt")),
0416       dRZHighEtThres_(pset.getParameter<std::vector<double> >("dRZMaxHighEtThres")),
0417       dRZLowEtGrad_(pset.getParameter<std::vector<double> >("dRZMaxLowEtGrad")),
0418       etaBins_(pset.getParameter<std::vector<double> >("etaBins")) {
0419   auto binSizeCheck = [](size_t sizeEtaBins, const std::vector<double>& vec, const std::string& name) {
0420     if (vec.size() != sizeEtaBins + 1) {
0421       throw cms::Exception("InvalidConfig")
0422           << " when constructing TrajSeedMatcher::MatchingCutsV2 " << name << " has " << vec.size()
0423           << " bins, it should be equal to #bins of etaBins+1" << sizeEtaBins + 1;
0424     }
0425   };
0426   binSizeCheck(etaBins_.size(), dPhiHighEt_, "dPhiMaxHighEt");
0427   binSizeCheck(etaBins_.size(), dPhiHighEtThres_, "dPhiMaxHighEtThres");
0428   binSizeCheck(etaBins_.size(), dPhiLowEtGrad_, "dPhiMaxLowEtGrad");
0429   binSizeCheck(etaBins_.size(), dRZHighEt_, "dRZMaxHighEt");
0430   binSizeCheck(etaBins_.size(), dRZHighEtThres_, "dRZMaxHighEtThres");
0431   binSizeCheck(etaBins_.size(), dRZLowEtGrad_, "dRZMaxLowEtGrad");
0432 }
0433 
0434 bool TrajSeedMatcher::MatchingCutsV2::operator()(const TrajSeedMatcher::SCHitMatch& scHitMatch) const {
0435   size_t binNr = getBinNr(scHitMatch.eta);
0436   float dPhiMax = getCutValue(scHitMatch.et, dPhiHighEt_[binNr], dPhiHighEtThres_[binNr], dPhiLowEtGrad_[binNr]);
0437   if (dPhiMax >= 0 && std::abs(scHitMatch.dPhi) > dPhiMax)
0438     return false;
0439   float dRZMax = getCutValue(scHitMatch.et, dRZHighEt_[binNr], dRZHighEtThres_[binNr], dRZLowEtGrad_[binNr]);
0440   if (dRZMax >= 0 && std::abs(scHitMatch.dRZ) > dRZMax)
0441     return false;
0442 
0443   return true;
0444 }
0445 
0446 //eta bins is exactly 1 smaller than the vectors which will be accessed by this bin nr
0447 size_t TrajSeedMatcher::MatchingCutsV2::getBinNr(float eta) const {
0448   const float absEta = std::abs(eta);
0449   for (size_t etaNr = 0; etaNr < etaBins_.size(); etaNr++) {
0450     if (absEta < etaBins_[etaNr])
0451       return etaNr;
0452   }
0453   return etaBins_.size();
0454 }