File indexing completed on 2021-10-01 22:40:59
0001 #include "RecoPixelVertexing/PixelTriplets/plugins/PixelTripletLargeTipGenerator.h"
0002 #include "RecoTracker/TkHitPairs/interface/HitPairGeneratorFromLayerPair.h"
0003 #include "FWCore/Framework/interface/ConsumesCollector.h"
0004
0005 #include "RecoPixelVertexing/PixelTriplets/interface/ThirdHitPredictionFromCircle.h"
0006 #include "RecoPixelVertexing/PixelTriplets/interface/ThirdHitRZPrediction.h"
0007 #include "RecoTracker/TkMSParametrization/interface/PixelRecoUtilities.h"
0008 #include "FWCore/Framework/interface/ConsumesCollector.h"
0009 #include "FWCore/Framework/interface/EventSetup.h"
0010 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0011 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0012
0013 #include "RecoPixelVertexing/PixelTriplets/plugins/ThirdHitCorrection.h"
0014 #include "RecoTracker/TkHitPairs/interface/RecHitsSortedInPhi.h"
0015
0016 #include "MatchedHitRZCorrectionFromBending.h"
0017 #include "CommonTools/RecoAlgos/interface/KDTreeLinkerAlgo.h"
0018
0019 #include <algorithm>
0020 #include <iostream>
0021 #include <vector>
0022 #include <cmath>
0023 #include <map>
0024
0025 #include "DataFormats/Math/interface/normalizedPhi.h"
0026
0027 #include "CommonTools/Utils/interface/DynArray.h"
0028
0029 using namespace std;
0030
0031 using Range = PixelRecoRange<float>;
0032 using HelixRZ = ThirdHitPredictionFromCircle::HelixRZ;
0033
0034 namespace {
0035 struct LayerRZPredictions {
0036 ThirdHitRZPrediction<PixelRecoLineRZ> line;
0037 ThirdHitRZPrediction<HelixRZ> helix1, helix2;
0038 MatchedHitRZCorrectionFromBending rzPositionFixup;
0039 ThirdHitCorrection correction;
0040 };
0041 }
0042
0043 constexpr double nSigmaRZ = 3.4641016151377544;
0044 constexpr float nSigmaPhi = 3.;
0045 constexpr float fnSigmaRZ = nSigmaRZ;
0046
0047 PixelTripletLargeTipGenerator::PixelTripletLargeTipGenerator(const edm::ParameterSet& cfg, edm::ConsumesCollector& iC)
0048 : HitTripletGeneratorFromPairAndLayers(cfg),
0049 useFixedPreFiltering(cfg.getParameter<bool>("useFixedPreFiltering")),
0050 extraHitRZtolerance(cfg.getParameter<double>("extraHitRZtolerance")),
0051 extraHitRPhitolerance(cfg.getParameter<double>("extraHitRPhitolerance")),
0052 useMScat(cfg.getParameter<bool>("useMultScattering")),
0053 useBend(cfg.getParameter<bool>("useBending")),
0054 dphi(useFixedPreFiltering ? cfg.getParameter<double>("phiPreFiltering") : 0),
0055 trackerTopologyESToken_(iC.esConsumes()),
0056 fieldESToken_(iC.esConsumes()) {
0057 if (useMScat) {
0058 msmakerESToken_ = iC.esConsumes();
0059 }
0060 }
0061
0062 PixelTripletLargeTipGenerator::~PixelTripletLargeTipGenerator() {}
0063
0064 void PixelTripletLargeTipGenerator::fillDescriptions(edm::ParameterSetDescription& desc) {
0065 HitTripletGeneratorFromPairAndLayers::fillDescriptions(desc);
0066
0067
0068 desc.add<double>("extraHitRPhitolerance", 0);
0069 desc.add<double>("extraHitRZtolerance", 0);
0070 desc.add<bool>("useMultScattering", true);
0071 desc.add<bool>("useBending", true);
0072 desc.add<bool>("useFixedPreFiltering", false);
0073 desc.add<double>("phiPreFiltering", 0.3);
0074 }
0075
0076 namespace {
0077 inline bool intersect(Range& range, const Range& second) {
0078 if ((range.min() > second.max()) | (range.max() < second.min()))
0079 return false;
0080 if (range.first < second.min())
0081 range.first = second.min();
0082 if (range.second > second.max())
0083 range.second = second.max();
0084 return range.first < range.second;
0085 }
0086 }
0087
0088 void PixelTripletLargeTipGenerator::hitTriplets(const TrackingRegion& region,
0089 OrderedHitTriplets& result,
0090 const edm::Event& ev,
0091 const edm::EventSetup& es,
0092 const SeedingLayerSetsHits::SeedingLayerSet& pairLayers,
0093 const std::vector<SeedingLayerSetsHits::SeedingLayer>& thirdLayers) {
0094 auto const& doublets = thePairGenerator->doublets(region, ev, es, pairLayers);
0095
0096 if (doublets.empty())
0097 return;
0098
0099 assert(theLayerCache);
0100 hitTriplets(region, result, ev, es, doublets, thirdLayers, nullptr, *theLayerCache);
0101 }
0102
0103 void PixelTripletLargeTipGenerator::hitTriplets(const TrackingRegion& region,
0104 OrderedHitTriplets& result,
0105 const edm::Event& ev,
0106 const edm::EventSetup& es,
0107 const HitDoublets& doublets,
0108 const std::vector<SeedingLayerSetsHits::SeedingLayer>& thirdLayers,
0109 std::vector<int>* tripletLastLayerIndex,
0110 LayerCacheType& layerCache) {
0111 int size = thirdLayers.size();
0112 const RecHitsSortedInPhi* thirdHitMap[size];
0113 vector<const DetLayer*> thirdLayerDetLayer(size, nullptr);
0114 for (int il = 0; il < size; ++il) {
0115 thirdHitMap[il] = &layerCache(thirdLayers[il], region);
0116 thirdLayerDetLayer[il] = thirdLayers[il].detLayer();
0117 }
0118 hitTriplets(region, result, es, doublets, thirdHitMap, thirdLayerDetLayer, size, tripletLastLayerIndex);
0119 }
0120
0121 void PixelTripletLargeTipGenerator::hitTriplets(const TrackingRegion& region,
0122 OrderedHitTriplets& result,
0123 const edm::EventSetup& es,
0124 const HitDoublets& doublets,
0125 const RecHitsSortedInPhi** thirdHitMap,
0126 const std::vector<const DetLayer*>& thirdLayerDetLayer,
0127 const int nThirdLayers) {
0128 hitTriplets(region, result, es, doublets, thirdHitMap, thirdLayerDetLayer, nThirdLayers, nullptr);
0129 }
0130
0131 void PixelTripletLargeTipGenerator::hitTriplets(const TrackingRegion& region,
0132 OrderedHitTriplets& result,
0133 const edm::EventSetup& es,
0134 const HitDoublets& doublets,
0135 const RecHitsSortedInPhi** thirdHitMap,
0136 const std::vector<const DetLayer*>& thirdLayerDetLayer,
0137 const int nThirdLayers,
0138 std::vector<int>* tripletLastLayerIndex) {
0139 const TrackerTopology* tTopo = &es.getData(trackerTopologyESToken_);
0140 const auto& field = es.getData(fieldESToken_);
0141 const MultipleScatteringParametrisationMaker* msmaker = nullptr;
0142 if (useMScat) {
0143 msmaker = &es.getData(msmakerESToken_);
0144 }
0145
0146 auto outSeq = doublets.detLayer(HitDoublets::outer)->seqNum();
0147
0148 using NodeInfo = KDTreeNodeInfo<unsigned int, 2>;
0149 std::vector<NodeInfo> layerTree;
0150 std::vector<unsigned int> foundNodes;
0151 foundNodes.reserve(100);
0152
0153 declareDynArray(KDTreeLinkerAlgo<unsigned int>, nThirdLayers, hitTree);
0154 declareDynArray(LayerRZPredictions, nThirdLayers, mapPred);
0155
0156 float rzError[nThirdLayers];
0157
0158 const float maxDelphi = region.ptMin() < 0.3f ? float(M_PI) / 4.f : float(M_PI) / 8.f;
0159 const float maxphi = M_PI + maxDelphi, minphi = -maxphi;
0160 const float safePhi = M_PI - maxDelphi;
0161
0162 for (int il = 0; il < nThirdLayers; il++) {
0163 auto const& hits = *thirdHitMap[il];
0164
0165 const DetLayer* layer = thirdLayerDetLayer[il];
0166 LayerRZPredictions& predRZ = mapPred[il];
0167 predRZ.line.initLayer(layer);
0168 predRZ.helix1.initLayer(layer);
0169 predRZ.helix2.initLayer(layer);
0170 predRZ.line.initTolerance(extraHitRZtolerance);
0171 predRZ.helix1.initTolerance(extraHitRZtolerance);
0172 predRZ.helix2.initTolerance(extraHitRZtolerance);
0173 predRZ.rzPositionFixup = MatchedHitRZCorrectionFromBending(layer, tTopo);
0174 predRZ.correction.init(region.ptMin(),
0175 *doublets.detLayer(HitDoublets::inner),
0176 *doublets.detLayer(HitDoublets::outer),
0177 *thirdLayerDetLayer[il],
0178 useMScat,
0179 msmaker,
0180 false,
0181 nullptr);
0182
0183 layerTree.clear();
0184 float minv = 999999.0;
0185 float maxv = -999999.0;
0186 float maxErr = 0.0f;
0187 for (unsigned int i = 0; i != hits.size(); ++i) {
0188 auto angle = hits.phi(i);
0189 auto v = hits.gv(i);
0190
0191 minv = std::min(minv, v);
0192 maxv = std::max(maxv, v);
0193 float myerr = hits.dv[i];
0194 maxErr = std::max(maxErr, myerr);
0195 layerTree.emplace_back(i, angle, v);
0196
0197 if (angle > safePhi)
0198 layerTree.emplace_back(i, angle - Geom::ftwoPi(), v);
0199 else if (angle < -safePhi)
0200 layerTree.emplace_back(i, angle + Geom::ftwoPi(), v);
0201 }
0202 KDTreeBox phiZ(minphi, maxphi, minv - 0.01f, maxv + 0.01f);
0203
0204 hitTree[il].build(layerTree, phiZ);
0205 rzError[il] = maxErr;
0206 }
0207
0208 double curv = PixelRecoUtilities::curvature(1. / region.ptMin(), field);
0209
0210 for (std::size_t ip = 0; ip != doublets.size(); ip++) {
0211 auto xi = doublets.x(ip, HitDoublets::inner);
0212 auto yi = doublets.y(ip, HitDoublets::inner);
0213 auto zi = doublets.z(ip, HitDoublets::inner);
0214
0215 auto xo = doublets.x(ip, HitDoublets::outer);
0216 auto yo = doublets.y(ip, HitDoublets::outer);
0217 auto zo = doublets.z(ip, HitDoublets::outer);
0218
0219 GlobalPoint gp1(xi, yi, zi);
0220 GlobalPoint gp2(xo, yo, zo);
0221
0222 auto toPos = std::signbit(zo - zi);
0223
0224 PixelRecoLineRZ line(gp1, gp2);
0225 PixelRecoPointRZ point2(gp2.perp(), zo);
0226 ThirdHitPredictionFromCircle predictionRPhi(gp1, gp2, extraHitRPhitolerance);
0227
0228 Range generalCurvature = predictionRPhi.curvature(region.originRBound());
0229 if (!intersect(generalCurvature, Range(-curv, curv)))
0230 continue;
0231
0232 for (int il = 0; il < nThirdLayers; il++) {
0233 if (hitTree[il].empty())
0234 continue;
0235 const DetLayer* layer = thirdLayerDetLayer[il];
0236 bool barrelLayer = layer->isBarrel();
0237
0238 if ((!barrelLayer) & (toPos != std::signbit(layer->position().z())))
0239 continue;
0240
0241 Range curvature = generalCurvature;
0242
0243 LayerRZPredictions& predRZ = mapPred[il];
0244 predRZ.line.initPropagator(&line);
0245
0246 auto& correction = predRZ.correction;
0247 correction.init(line, point2, outSeq);
0248
0249 Range rzRange;
0250 if (useBend) {
0251
0252
0253
0254
0255
0256
0257
0258
0259
0260
0261
0262
0263
0264
0265
0266
0267
0268
0269 if (!barrelLayer) {
0270 Range z3s = predRZ.line.detRange();
0271 double z3 = z3s.first < 0 ? std::max(z3s.first, z3s.second) : std::min(z3s.first, z3s.second);
0272 double maxCurvature = HelixRZ::maxCurvature(&predictionRPhi, gp1.z(), gp2.z(), z3);
0273 if (!intersect(curvature, Range(-maxCurvature, maxCurvature)))
0274 continue;
0275 }
0276
0277 HelixRZ helix1(&predictionRPhi, gp1.z(), gp2.z(), curvature.first);
0278 HelixRZ helix2(&predictionRPhi, gp1.z(), gp2.z(), curvature.second);
0279
0280 predRZ.helix1.initPropagator(&helix1);
0281 predRZ.helix2.initPropagator(&helix2);
0282
0283 Range rzRanges[2] = {predRZ.helix1(), predRZ.helix2()};
0284 predRZ.helix1.initPropagator(nullptr);
0285 predRZ.helix2.initPropagator(nullptr);
0286
0287 rzRange.first = std::min(rzRanges[0].first, rzRanges[1].first);
0288 rzRange.second = std::max(rzRanges[0].second, rzRanges[1].second);
0289
0290
0291
0292 if (curvature.first * curvature.second < 0.0) {
0293 Range rzLineRange = predRZ.line();
0294 rzRange.first = std::min(rzRange.first, rzLineRange.first);
0295 rzRange.second = std::max(rzRange.second, rzLineRange.second);
0296 }
0297 } else {
0298 rzRange = predRZ.line();
0299 }
0300
0301 if (rzRange.first >= rzRange.second)
0302 continue;
0303
0304 correction.correctRZRange(rzRange);
0305
0306 Range phiRange;
0307 if (useFixedPreFiltering) {
0308 float phi0 = doublets.phi(ip, HitDoublets::outer);
0309 phiRange = Range(phi0 - dphi, phi0 + dphi);
0310 } else {
0311 Range radius;
0312
0313 if (barrelLayer) {
0314 radius = predRZ.line.detRange();
0315 if (!intersect(rzRange, predRZ.line.detSize()))
0316 continue;
0317 } else {
0318 radius = rzRange;
0319 if (!intersect(radius, predRZ.line.detSize()))
0320 continue;
0321 }
0322
0323
0324 if ((curvature.first < 0.0f) & (curvature.second < 0.0f)) {
0325 radius.swap();
0326 } else if ((curvature.first >= 0.0f) & (curvature.second >= 0.0f)) {
0327 ;
0328 } else {
0329 radius.first = radius.second;
0330 }
0331 auto phi12 = predictionRPhi.phi(curvature.first, radius.first);
0332 auto phi22 = predictionRPhi.phi(curvature.second, radius.second);
0333 phi12 = normalizedPhi(phi12);
0334 phi22 = proxim(phi22, phi12);
0335 phiRange = Range(phi12, phi22);
0336 phiRange.sort();
0337 auto rmean = radius.mean();
0338 phiRange.first *= rmean;
0339 phiRange.second *= rmean;
0340 correction.correctRPhiRange(phiRange);
0341 phiRange.first /= rmean;
0342 phiRange.second /= rmean;
0343 }
0344
0345 foundNodes.clear();
0346 float prmin = phiRange.min(), prmax = phiRange.max();
0347
0348 if (prmax - prmin > maxDelphi) {
0349 auto prm = phiRange.mean();
0350 prmin = prm - 0.5f * maxDelphi;
0351 prmax = prm + 0.5f * maxDelphi;
0352 }
0353
0354 if (barrelLayer) {
0355 Range regMax = predRZ.line.detRange();
0356 Range regMin = predRZ.line(regMax.min());
0357 regMax = predRZ.line(regMax.max());
0358 correction.correctRZRange(regMin);
0359 correction.correctRZRange(regMax);
0360 if (regMax.min() < regMin.min()) {
0361 std::swap(regMax, regMin);
0362 }
0363 KDTreeBox phiZ(prmin, prmax, regMin.min() - fnSigmaRZ * rzError[il], regMax.max() + fnSigmaRZ * rzError[il]);
0364 hitTree[il].search(phiZ, foundNodes);
0365 } else {
0366 KDTreeBox phiZ(prmin, prmax, rzRange.min() - fnSigmaRZ * rzError[il], rzRange.max() + fnSigmaRZ * rzError[il]);
0367 hitTree[il].search(phiZ, foundNodes);
0368 }
0369
0370 MatchedHitRZCorrectionFromBending l2rzFixup(doublets.hit(ip, HitDoublets::outer)->det()->geographicalId(), tTopo);
0371 MatchedHitRZCorrectionFromBending l3rzFixup = predRZ.rzPositionFixup;
0372
0373 auto const& hits = *thirdHitMap[il];
0374 for (auto KDdata : foundNodes) {
0375 GlobalPoint p3 = hits.gp(KDdata);
0376 double p3_r = p3.perp();
0377 double p3_z = p3.z();
0378 float p3_phi = hits.phi(KDdata);
0379
0380 Range rangeRPhi = predictionRPhi(curvature, p3_r);
0381 correction.correctRPhiRange(rangeRPhi);
0382
0383 float ir = 1.f / p3_r;
0384
0385 constexpr float maxPhiErr = 0.5 * M_PI;
0386 float phiErr = nSigmaPhi * hits.drphi[KDdata] * ir;
0387 phiErr = std::min(maxPhiErr, phiErr);
0388 if (!checkPhiInRange(p3_phi, rangeRPhi.first * ir - phiErr, rangeRPhi.second * ir + phiErr, maxPhiErr))
0389 continue;
0390
0391 Basic2DVector<double> thc(p3.x(), p3.y());
0392
0393 auto curv_ = predictionRPhi.curvature(thc);
0394 double p2_r = point2.r();
0395 double p2_z = point2.z();
0396
0397 l2rzFixup(predictionRPhi, curv_, *doublets.hit(ip, HitDoublets::outer), p2_r, p2_z, tTopo);
0398 l3rzFixup(predictionRPhi, curv_, *hits.theHits[KDdata].hit(), p3_r, p3_z, tTopo);
0399
0400 Range rangeRZ;
0401 if (useBend) {
0402 HelixRZ updatedHelix(&predictionRPhi, gp1.z(), p2_z, curv_);
0403 rangeRZ = predRZ.helix1(barrelLayer ? p3_r : p3_z, updatedHelix);
0404 } else {
0405 float tIP = predictionRPhi.transverseIP(thc);
0406 PixelRecoPointRZ updatedPoint2(p2_r, p2_z);
0407 PixelRecoLineRZ updatedLine(line.origin(), point2, tIP);
0408 rangeRZ = predRZ.line(barrelLayer ? p3_r : p3_z, line);
0409 }
0410 correction.correctRZRange(rangeRZ);
0411
0412 double err = nSigmaRZ * hits.dv[KDdata];
0413
0414 rangeRZ.first -= err, rangeRZ.second += err;
0415
0416 if (!rangeRZ.inside(barrelLayer ? p3_z : p3_r))
0417 continue;
0418
0419 if (theMaxElement != 0 && result.size() >= theMaxElement) {
0420 result.clear();
0421 if (tripletLastLayerIndex)
0422 tripletLastLayerIndex->clear();
0423 edm::LogError("TooManyTriplets") << " number of triples exceed maximum. no triplets produced.";
0424 return;
0425 }
0426 result.emplace_back(
0427 doublets.hit(ip, HitDoublets::inner), doublets.hit(ip, HitDoublets::outer), hits.theHits[KDdata].hit());
0428
0429 if (tripletLastLayerIndex)
0430 tripletLastLayerIndex->push_back(il);
0431 }
0432 }
0433 }
0434
0435 }