Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:24:48 

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   for (auto layer : inLayers) {
0327     if (GeomDetEnumerators::isTrackerPixel(layer->subDetector())) {
0328       if (layerHasValidHits(*layer, hitTrajState, backwardPropagator_))
0329         nrValidLayers++;
0330     }
0331   }
0332   for (auto layer : outLayers) {
0333     if (GeomDetEnumerators::isTrackerPixel(layer->subDetector())) {
0334       if (layerHasValidHits(*layer, hitTrajState, forwardPropagator_))
0335         nrValidLayers++;
0336     }
0337   }
0338   return nrValidLayers;
0339 }
0340 
0341 bool TrajSeedMatcher::layerHasValidHits(const DetLayer& layer,
0342                                         const TrajectoryStateOnSurface& hitSurState,
0343                                         const Propagator& propToLayerFromState) const {
0344   //FIXME: do not hardcode with weird magic numbers stolen from ancient tracking code
0345   //its taken from https://cmssdt.cern.ch/dxr/CMSSW/source/RecoTracker/TrackProducer/interface/TrackProducerBase.icc#165
0346   //which inspires this code
0347   Chi2MeasurementEstimator estimator(30., -3.0, 0.5, 2.0, 0.5, 1.e12);  // same as defauts....
0348 
0349   const std::vector<GeometricSearchDet::DetWithState>& detWithState =
0350       layer.compatibleDets(hitSurState, propToLayerFromState, estimator);
0351   if (detWithState.empty())
0352     return false;
0353   else {
0354     DetId id = detWithState.front().first->geographicalId();
0355     MeasurementDetWithData measDet = measTkEvt_.idToDet(id);
0356     //Below, measDet.hasBadComponents handles the check that a Pixel module has or not errors like FED25
0357     if (measDet.isActive() && !measDet.hasBadComponents(detWithState.front().second))
0358       return true;
0359     else
0360       return false;
0361   }
0362 }
0363 
0364 size_t TrajSeedMatcher::getNrHitsRequired(const int nrValidLayers) const {
0365   for (size_t binNr = 0; binNr < cfg_.minNrHitsValidLayerBins.size(); binNr++) {
0366     if (nrValidLayers < cfg_.minNrHitsValidLayerBins[binNr])
0367       return cfg_.minNrHits[binNr];
0368   }
0369   return cfg_.minNrHits.back();
0370 }
0371 
0372 TrajSeedMatcher::MatchingCutsV1::MatchingCutsV1(const edm::ParameterSet& pset)
0373     : dPhiMax_(pset.getParameter<double>("dPhiMax")),
0374       dRZMax_(pset.getParameter<double>("dRZMax")),
0375       dRZMaxLowEtThres_(pset.getParameter<double>("dRZMaxLowEtThres")),
0376       dRZMaxLowEtEtaBins_(pset.getParameter<std::vector<double> >("dRZMaxLowEtEtaBins")),
0377       dRZMaxLowEt_(pset.getParameter<std::vector<double> >("dRZMaxLowEt")) {
0378   if (dRZMaxLowEtEtaBins_.size() + 1 != dRZMaxLowEt_.size()) {
0379     throw cms::Exception("InvalidConfig") << " dRZMaxLowEtEtaBins should be 1 less than dRZMaxLowEt when its "
0380                                           << dRZMaxLowEtEtaBins_.size() << " vs " << dRZMaxLowEt_.size();
0381   }
0382 }
0383 
0384 bool TrajSeedMatcher::MatchingCutsV1::operator()(const TrajSeedMatcher::SCHitMatch& scHitMatch) const {
0385   if (dPhiMax_ >= 0 && std::abs(scHitMatch.dPhi) > dPhiMax_)
0386     return false;
0387 
0388   const float dRZMax = getDRZCutValue(scHitMatch.et, scHitMatch.eta);
0389   if (dRZMax_ >= 0 && std::abs(scHitMatch.dRZ) > dRZMax)
0390     return false;
0391 
0392   return true;
0393 }
0394 
0395 float TrajSeedMatcher::MatchingCutsV1::getDRZCutValue(const float scEt, const float scEta) const {
0396   if (scEt >= dRZMaxLowEtThres_)
0397     return dRZMax_;
0398   else {
0399     const float absEta = std::abs(scEta);
0400     for (size_t etaNr = 0; etaNr < dRZMaxLowEtEtaBins_.size(); etaNr++) {
0401       if (absEta < dRZMaxLowEtEtaBins_[etaNr])
0402         return dRZMaxLowEt_[etaNr];
0403     }
0404     return dRZMaxLowEt_.back();
0405   }
0406 }
0407 
0408 TrajSeedMatcher::MatchingCutsV2::MatchingCutsV2(const edm::ParameterSet& pset)
0409     : dPhiHighEt_(pset.getParameter<std::vector<double> >("dPhiMaxHighEt")),
0410       dPhiHighEtThres_(pset.getParameter<std::vector<double> >("dPhiMaxHighEtThres")),
0411       dPhiLowEtGrad_(pset.getParameter<std::vector<double> >("dPhiMaxLowEtGrad")),
0412       dRZHighEt_(pset.getParameter<std::vector<double> >("dRZMaxHighEt")),
0413       dRZHighEtThres_(pset.getParameter<std::vector<double> >("dRZMaxHighEtThres")),
0414       dRZLowEtGrad_(pset.getParameter<std::vector<double> >("dRZMaxLowEtGrad")),
0415       etaBins_(pset.getParameter<std::vector<double> >("etaBins")) {
0416   auto binSizeCheck = [](size_t sizeEtaBins, const std::vector<double>& vec, const std::string& name) {
0417     if (vec.size() != sizeEtaBins + 1) {
0418       throw cms::Exception("InvalidConfig")
0419           << " when constructing TrajSeedMatcher::MatchingCutsV2 " << name << " has " << vec.size()
0420           << " bins, it should be equal to #bins of etaBins+1" << sizeEtaBins + 1;
0421     }
0422   };
0423   binSizeCheck(etaBins_.size(), dPhiHighEt_, "dPhiMaxHighEt");
0424   binSizeCheck(etaBins_.size(), dPhiHighEtThres_, "dPhiMaxHighEtThres");
0425   binSizeCheck(etaBins_.size(), dPhiLowEtGrad_, "dPhiMaxLowEtGrad");
0426   binSizeCheck(etaBins_.size(), dRZHighEt_, "dRZMaxHighEt");
0427   binSizeCheck(etaBins_.size(), dRZHighEtThres_, "dRZMaxHighEtThres");
0428   binSizeCheck(etaBins_.size(), dRZLowEtGrad_, "dRZMaxLowEtGrad");
0429 }
0430 
0431 bool TrajSeedMatcher::MatchingCutsV2::operator()(const TrajSeedMatcher::SCHitMatch& scHitMatch) const {
0432   size_t binNr = getBinNr(scHitMatch.eta);
0433   float dPhiMax = getCutValue(scHitMatch.et, dPhiHighEt_[binNr], dPhiHighEtThres_[binNr], dPhiLowEtGrad_[binNr]);
0434   if (dPhiMax >= 0 && std::abs(scHitMatch.dPhi) > dPhiMax)
0435     return false;
0436   float dRZMax = getCutValue(scHitMatch.et, dRZHighEt_[binNr], dRZHighEtThres_[binNr], dRZLowEtGrad_[binNr]);
0437   if (dRZMax >= 0 && std::abs(scHitMatch.dRZ) > dRZMax)
0438     return false;
0439 
0440   return true;
0441 }
0442 
0443 //eta bins is exactly 1 smaller than the vectors which will be accessed by this bin nr
0444 size_t TrajSeedMatcher::MatchingCutsV2::getBinNr(float eta) const {
0445   const float absEta = std::abs(eta);
0446   for (size_t etaNr = 0; etaNr < etaBins_.size(); etaNr++) {
0447     if (absEta < etaBins_[etaNr])
0448       return etaNr;
0449   }
0450   return etaBins_.size();
0451 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!