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