File indexing completed on 2024-10-08 23:10:01
0001 #include "RecoTracker/PixelSeeding/plugins/PixelTripletHLTGenerator.h"
0002 #include "RecoTracker/TkHitPairs/interface/HitPairGeneratorFromLayerPair.h"
0003 #include "FWCore/Framework/interface/ConsumesCollector.h"
0004 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0005
0006 #include "ThirdHitPredictionFromInvParabola.h"
0007 #include "RecoTracker/PixelSeeding/interface/ThirdHitRZPrediction.h"
0008 #include "RecoTracker/TkMSParametrization/interface/PixelRecoUtilities.h"
0009
0010 #include "ThirdHitCorrection.h"
0011 #include "RecoTracker/TkHitPairs/interface/RecHitsSortedInPhi.h"
0012 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0013
0014 #include "RecoTracker/TkSeedingLayers/interface/SeedComparitorFactory.h"
0015 #include "RecoTracker/TkSeedingLayers/interface/SeedComparitor.h"
0016
0017 #include "DataFormats/GeometryVector/interface/Pi.h"
0018 #include "CommonTools/RecoAlgos/interface/KDTreeLinkerAlgo.h"
0019
0020 #include "CommonTools/Utils/interface/DynArray.h"
0021
0022 #include "DataFormats/Math/interface/normalizedPhi.h"
0023
0024 #include <cstdio>
0025 #include <iostream>
0026
0027 using pixelrecoutilities::LongitudinalBendingCorrection;
0028 using Range = PixelRecoRange<float>;
0029
0030 using namespace std;
0031
0032 PixelTripletHLTGenerator::PixelTripletHLTGenerator(const edm::ParameterSet& cfg, edm::ConsumesCollector& iC)
0033 : HitTripletGeneratorFromPairAndLayers(cfg),
0034 fieldToken_(iC.esConsumes()),
0035 useFixedPreFiltering(cfg.getParameter<bool>("useFixedPreFiltering")),
0036 extraHitRZtolerance(cfg.getParameter<double>("extraHitRZtolerance")),
0037 extraHitRPhitolerance(cfg.getParameter<double>("extraHitRPhitolerance")),
0038 useMScat(cfg.getParameter<bool>("useMultScattering")),
0039 useBend(cfg.getParameter<bool>("useBending")),
0040 dphi(useFixedPreFiltering ? cfg.getParameter<double>("phiPreFiltering") : 0) {
0041 edm::ParameterSet comparitorPSet = cfg.getParameter<edm::ParameterSet>("SeedComparitorPSet");
0042 std::string comparitorName = comparitorPSet.getParameter<std::string>("ComponentName");
0043 if (comparitorName != "none") {
0044 theComparitor = SeedComparitorFactory::get()->create(comparitorName, comparitorPSet, iC);
0045 }
0046 if (useMScat) {
0047 msmakerToken_ = iC.esConsumes();
0048 }
0049 }
0050
0051 PixelTripletHLTGenerator::~PixelTripletHLTGenerator() {}
0052
0053 void PixelTripletHLTGenerator::fillDescriptions(edm::ParameterSetDescription& desc) {
0054 HitTripletGeneratorFromPairAndLayers::fillDescriptions(desc);
0055 desc.add<double>("extraHitRPhitolerance", 0.032);
0056 desc.add<double>("extraHitRZtolerance", 0.037);
0057 desc.add<bool>("useMultScattering", true);
0058 desc.add<bool>("useBending", true);
0059 desc.add<bool>("useFixedPreFiltering", false);
0060 desc.add<double>("phiPreFiltering", 0.3);
0061
0062 edm::ParameterSetDescription descComparitor;
0063 descComparitor.add<std::string>("ComponentName", "none");
0064 descComparitor.setAllowAnything();
0065 desc.add<edm::ParameterSetDescription>("SeedComparitorPSet", descComparitor);
0066 }
0067
0068 void PixelTripletHLTGenerator::hitTriplets(const TrackingRegion& region,
0069 OrderedHitTriplets& result,
0070 const edm::Event& ev,
0071 const edm::EventSetup& es,
0072 const SeedingLayerSetsHits::SeedingLayerSet& pairLayers,
0073 const std::vector<SeedingLayerSetsHits::SeedingLayer>& thirdLayers) {
0074 auto const& doublets = thePairGenerator->doublets(region, ev, es, pairLayers);
0075 if (not doublets or doublets->empty())
0076 return;
0077
0078 assert(theLayerCache);
0079 hitTriplets(region, result, ev, es, *doublets, thirdLayers, nullptr, *theLayerCache);
0080 }
0081
0082 void PixelTripletHLTGenerator::hitTriplets(const TrackingRegion& region,
0083 OrderedHitTriplets& result,
0084 const edm::Event& ev,
0085 const edm::EventSetup& es,
0086 const HitDoublets& doublets,
0087 const std::vector<SeedingLayerSetsHits::SeedingLayer>& thirdLayers,
0088 std::vector<int>* tripletLastLayerIndex,
0089 LayerCacheType& layerCache) {
0090 if (theComparitor)
0091 theComparitor->init(ev, es);
0092
0093 int size = thirdLayers.size();
0094 const RecHitsSortedInPhi* thirdHitMap[size];
0095 vector<const DetLayer*> thirdLayerDetLayer(size, nullptr);
0096 for (int il = 0; il < size; ++il) {
0097 thirdHitMap[il] = &layerCache(thirdLayers[il], region);
0098 thirdLayerDetLayer[il] = thirdLayers[il].detLayer();
0099 }
0100 hitTriplets(region, result, es, doublets, thirdHitMap, thirdLayerDetLayer, size, tripletLastLayerIndex);
0101 }
0102
0103 void PixelTripletHLTGenerator::hitTriplets(const TrackingRegion& region,
0104 OrderedHitTriplets& result,
0105 const edm::EventSetup& es,
0106 const HitDoublets& doublets,
0107 const RecHitsSortedInPhi** thirdHitMap,
0108 const std::vector<const DetLayer*>& thirdLayerDetLayer,
0109 const int nThirdLayers) {
0110 hitTriplets(region, result, es, doublets, thirdHitMap, thirdLayerDetLayer, nThirdLayers, nullptr);
0111 }
0112
0113 void PixelTripletHLTGenerator::hitTriplets(const TrackingRegion& region,
0114 OrderedHitTriplets& result,
0115 const edm::EventSetup& es,
0116 const HitDoublets& doublets,
0117 const RecHitsSortedInPhi** thirdHitMap,
0118 const std::vector<const DetLayer*>& thirdLayerDetLayer,
0119 const int nThirdLayers,
0120 std::vector<int>* tripletLastLayerIndex) {
0121 auto outSeq = doublets.detLayer(HitDoublets::outer)->seqNum();
0122
0123 float regOffset = region.origin().perp();
0124
0125 declareDynArray(ThirdHitRZPrediction<PixelRecoLineRZ>, nThirdLayers, preds);
0126 declareDynArray(ThirdHitCorrection, nThirdLayers, corrections);
0127
0128 typedef RecHitsSortedInPhi::Hit Hit;
0129
0130 using NodeInfo = KDTreeNodeInfo<unsigned int, 2>;
0131 std::vector<NodeInfo> layerTree;
0132 std::vector<unsigned int> foundNodes;
0133 foundNodes.reserve(100);
0134
0135 declareDynArray(KDTreeLinkerAlgo<unsigned int>, nThirdLayers, hitTree);
0136 std::vector<float> rzError(nThirdLayers);
0137
0138 const float maxDelphi = region.ptMin() < 0.3f ? float(M_PI) / 4.f : float(M_PI) / 8.f;
0139 const float maxphi = M_PI + maxDelphi, minphi = -maxphi;
0140 const float safePhi = M_PI - maxDelphi;
0141
0142 const auto& field = es.getData(fieldToken_);
0143 const MultipleScatteringParametrisationMaker* msmaker = nullptr;
0144 if (useMScat) {
0145 msmaker = &es.getData(msmakerToken_);
0146 }
0147
0148
0149 for (int il = 0; il < nThirdLayers; ++il) {
0150 auto const& hits = *thirdHitMap[il];
0151 ThirdHitRZPrediction<PixelRecoLineRZ>& pred = preds[il];
0152 pred.initLayer(thirdLayerDetLayer[il]);
0153 pred.initTolerance(extraHitRZtolerance);
0154
0155 corrections[il].init(region.ptMin(),
0156 *doublets.detLayer(HitDoublets::inner),
0157 *doublets.detLayer(HitDoublets::outer),
0158 *thirdLayerDetLayer[il],
0159 useMScat,
0160 msmaker,
0161 useBend,
0162 &field);
0163
0164 layerTree.clear();
0165 float minv = 999999.0f, maxv = -minv;
0166 float maxErr = 0.0f;
0167 for (unsigned int i = 0; i != hits.size(); ++i) {
0168 auto angle = hits.phi(i);
0169 auto v = hits.gv(i);
0170
0171 minv = std::min(minv, v);
0172 maxv = std::max(maxv, v);
0173 float myerr = hits.dv[i];
0174 maxErr = std::max(maxErr, myerr);
0175 layerTree.emplace_back(i, angle, v);
0176
0177 if (angle > safePhi)
0178 layerTree.emplace_back(i, angle - Geom::ftwoPi(), v);
0179 else if (angle < -safePhi)
0180 layerTree.emplace_back(i, angle + Geom::ftwoPi(), v);
0181 }
0182 KDTreeBox phiZ(minphi, maxphi, minv - 0.01f, maxv + 0.01f);
0183
0184 hitTree[il].build(layerTree, phiZ);
0185 rzError[il] = maxErr;
0186
0187 }
0188
0189 float imppar = region.originRBound();
0190 float imppartmp = region.originRBound() + region.origin().perp();
0191 float curv = PixelRecoUtilities::curvature(1.f / region.ptMin(), field);
0192
0193 for (std::size_t ip = 0; ip != doublets.size(); ip++) {
0194 auto xi = doublets.x(ip, HitDoublets::inner);
0195 auto yi = doublets.y(ip, HitDoublets::inner);
0196 auto zi = doublets.z(ip, HitDoublets::inner);
0197 auto rvi = doublets.rv(ip, HitDoublets::inner);
0198 auto xo = doublets.x(ip, HitDoublets::outer);
0199 auto yo = doublets.y(ip, HitDoublets::outer);
0200 auto zo = doublets.z(ip, HitDoublets::outer);
0201 auto rvo = doublets.rv(ip, HitDoublets::outer);
0202
0203 auto toPos = std::signbit(zo - zi);
0204
0205 PixelRecoPointRZ point1(rvi, zi);
0206 PixelRecoPointRZ point2(rvo, zo);
0207 PixelRecoLineRZ line(point1, point2);
0208 ThirdHitPredictionFromInvParabola predictionRPhi(xi - region.origin().x(),
0209 yi - region.origin().y(),
0210 xo - region.origin().x(),
0211 yo - region.origin().y(),
0212 imppar,
0213 curv,
0214 extraHitRPhitolerance);
0215
0216 ThirdHitPredictionFromInvParabola predictionRPhitmp(xi, yi, xo, yo, imppartmp, curv, extraHitRPhitolerance);
0217
0218
0219
0220
0221
0222
0223 for (int il = 0; il != nThirdLayers; ++il) {
0224 const DetLayer* layer = thirdLayerDetLayer[il];
0225 auto barrelLayer = layer->isBarrel();
0226
0227 if ((!barrelLayer) & (toPos != std::signbit(layer->position().z())))
0228 continue;
0229
0230 if (hitTree[il].empty())
0231 continue;
0232
0233 auto const& hits = *thirdHitMap[il];
0234
0235 auto& correction = corrections[il];
0236
0237 correction.init(line, point2, outSeq);
0238
0239 auto& predictionRZ = preds[il];
0240
0241 predictionRZ.initPropagator(&line);
0242 Range rzRange = predictionRZ();
0243 correction.correctRZRange(rzRange);
0244
0245 Range phiRange;
0246 if (useFixedPreFiltering) {
0247 float phi0 = doublets.phi(ip, HitDoublets::outer);
0248 phiRange = Range(phi0 - dphi, phi0 + dphi);
0249 } else {
0250 Range radius;
0251 if (barrelLayer) {
0252 radius = predictionRZ.detRange();
0253 } else {
0254 radius =
0255 Range(max(rzRange.min(), predictionRZ.detSize().min()), min(rzRange.max(), predictionRZ.detSize().max()));
0256 }
0257 if (radius.empty())
0258 continue;
0259
0260
0261
0262 auto rPhi1 = predictionRPhitmp(radius.max());
0263 bool ok1 = !rPhi1.empty();
0264 if (ok1) {
0265 correction.correctRPhiRange(rPhi1);
0266 rPhi1.first /= radius.max();
0267 rPhi1.second /= radius.max();
0268 }
0269 auto rPhi2 = predictionRPhitmp(radius.min());
0270 bool ok2 = !rPhi2.empty();
0271 if (ok2) {
0272 correction.correctRPhiRange(rPhi2);
0273 rPhi2.first /= radius.min();
0274 rPhi2.second /= radius.min();
0275 }
0276
0277 if (ok1) {
0278 rPhi1.first = normalizedPhi(rPhi1.first);
0279 rPhi1.second = proxim(rPhi1.second, rPhi1.first);
0280 if (ok2) {
0281 rPhi2.first = proxim(rPhi2.first, rPhi1.first);
0282 rPhi2.second = proxim(rPhi2.second, rPhi1.first);
0283 phiRange = rPhi1.sum(rPhi2);
0284 } else
0285 phiRange = rPhi1;
0286 } else if (ok2) {
0287 rPhi2.first = normalizedPhi(rPhi2.first);
0288 rPhi2.second = proxim(rPhi2.second, rPhi2.first);
0289 phiRange = rPhi2;
0290 } else
0291 continue;
0292 }
0293
0294 constexpr float nSigmaRZ = 3.46410161514f;
0295 constexpr float nSigmaPhi = 3.f;
0296
0297 foundNodes.clear();
0298 float prmin = phiRange.min(), prmax = phiRange.max();
0299
0300 if (prmax - prmin > maxDelphi) {
0301 auto prm = phiRange.mean();
0302 prmin = prm - 0.5f * maxDelphi;
0303 prmax = prm + 0.5f * maxDelphi;
0304 }
0305
0306 if (barrelLayer) {
0307 Range regMax = predictionRZ.detRange();
0308 Range regMin = predictionRZ(regMax.min() - regOffset);
0309 regMax = predictionRZ(regMax.max() + regOffset);
0310 correction.correctRZRange(regMin);
0311 correction.correctRZRange(regMax);
0312 if (regMax.min() < regMin.min()) {
0313 swap(regMax, regMin);
0314 }
0315 KDTreeBox phiZ(prmin, prmax, regMin.min() - nSigmaRZ * rzError[il], regMax.max() + nSigmaRZ * rzError[il]);
0316 hitTree[il].search(phiZ, foundNodes);
0317 } else {
0318 KDTreeBox phiZ(prmin,
0319 prmax,
0320 rzRange.min() - regOffset - nSigmaRZ * rzError[il],
0321 rzRange.max() + regOffset + nSigmaRZ * rzError[il]);
0322 hitTree[il].search(phiZ, foundNodes);
0323 }
0324
0325
0326
0327
0328 for (auto KDdata : foundNodes) {
0329 if (theMaxElement != 0 && result.size() >= theMaxElement) {
0330 result.clear();
0331 if (tripletLastLayerIndex)
0332 tripletLastLayerIndex->clear();
0333 edm::LogError("TooManyTriplets") << " number of triples exceeds maximum. no triplets produced.";
0334 return;
0335 }
0336
0337 float p3_u = hits.u[KDdata];
0338 float p3_v = hits.v[KDdata];
0339 float p3_phi = hits.lphi[KDdata];
0340
0341
0342
0343
0344 Range allowed = predictionRZ(p3_u);
0345 correction.correctRZRange(allowed);
0346 float vErr = nSigmaRZ * hits.dv[KDdata];
0347 Range hitRange(p3_v - vErr, p3_v + vErr);
0348 Range crossingRange = allowed.intersection(hitRange);
0349 if (crossingRange.empty())
0350 continue;
0351
0352 float ir = 1.f / hits.rv(KDdata);
0353
0354 constexpr float maxPhiErr = 0.5 * M_PI;
0355 float phiErr = nSigmaPhi * hits.drphi[KDdata] * ir;
0356 phiErr = std::min(maxPhiErr, phiErr);
0357 bool nook = true;
0358 for (int icharge = -1; icharge <= 1; icharge += 2) {
0359 Range rangeRPhi = predictionRPhi(hits.rv(KDdata), icharge);
0360 if (rangeRPhi.first > rangeRPhi.second)
0361 continue;
0362 correction.correctRPhiRange(rangeRPhi);
0363 if (checkPhiInRange(p3_phi, rangeRPhi.first * ir - phiErr, rangeRPhi.second * ir + phiErr, maxPhiErr)) {
0364
0365 OrderedHitTriplet hittriplet(
0366 doublets.hit(ip, HitDoublets::inner), doublets.hit(ip, HitDoublets::outer), hits.theHits[KDdata].hit());
0367 if (!theComparitor || theComparitor->compatible(hittriplet)) {
0368 result.push_back(hittriplet);
0369
0370 if (tripletLastLayerIndex)
0371 tripletLastLayerIndex->push_back(il);
0372 } else {
0373 LogDebug("RejectedTriplet") << "rejected triplet from comparitor ";
0374 }
0375 nook = false;
0376 break;
0377 }
0378 }
0379 if (nook)
0380 LogDebug("RejectedTriplet") << "rejected triplet from second phicheck " << p3_phi;
0381 }
0382 }
0383 }
0384
0385 }