PixelSeeding/plugins/PixelTripletLargeTipGenerator.cc

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 }  // namespace
0042
0043 constexpr double nSigmaRZ = 3.4641016151377544;  // sqrt(12.)
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   // Defaults for the extraHit*tolerance are set to 0 here since that
0067   // was already the case in practice in all offline occurrances.
0068   desc.add<double>("extraHitRPhitolerance", 0);  // default in old python was 0.032
0069   desc.add<double>("extraHitRZtolerance", 0);    // default in old python was 0.037
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 // Disable bitwise-instead-of-logical warning, see discussion in
0077 // https://github.com/cms-sw/cmssw/issues/39105
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 }  // namespace
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;       // re-used throughout
0166   std::vector<unsigned int> foundNodes;  // re-used throughout
0167   foundNodes.reserve(100);
0168
0169   declareDynArray(KDTreeLinkerAlgo<unsigned int>, nThirdLayers, hitTree);
0170   declareDynArray(LayerRZPredictions, nThirdLayers, mapPred);
0171
0172   float rzError[nThirdLayers];  //save maximum errors
0173
0174   const float maxDelphi = region.ptMin() < 0.3f ? float(M_PI) / 4.f : float(M_PI) / 8.f;  // FIXME move to config??
0175   const float maxphi = M_PI + maxDelphi, minphi = -maxphi;  // increase to cater for any range
0176   const float safePhi = M_PI - maxDelphi;                   // sideband
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;  // Initialise to extreme values in case no hits
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       //use (phi,r) for endcaps rather than (phi,z)
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);  // save it
0212       // populate side-bands
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);  // declare our bounds
0219     //add fudge factors in case only one hit and also for floating-point inaccuracy
0220     hitTree[il].build(layerTree, phiZ);  // make KDtree
0221     rzError[il] = maxErr;                //save error
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     // auto rvi = doublets.rv(ip,HitDoublets::inner);
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     // auto rvo = doublets.rv(ip,HitDoublets::outer);
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;  // Don't bother if no hits
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         // For the barrel region:
0268         // swiping the helix passing through the two points across from
0269         // negative to positive bending, can give us a sort of U-shaped
0270         // projection onto the phi-z (barrel) or r-z plane (forward)
0271         // so we checking minimum/maximum of all three possible extrema
0272         //
0273         // For the endcap region:
0274         // Checking minimum/maximum radius of the helix projection
0275         // onto an endcap plane, here we have to guard against
0276         // looping tracks, when phi(delta z) gets out of control.
0277         // HelixRZ::rAtZ should not follow looping tracks, but clamp
0278         // to the minimum reachable r with the next-best lower |curvature|.
0279         // So same procedure as for the barrel region can be applied.
0280         //
0281         // In order to avoid looking for potential looping tracks at all
0282         // we also clamp the allowed curvature range for this layer,
0283         // and potentially fail the layer entirely
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         // if the allowed curvatures include a straight line,
0307         // this can give us another extremum for allowed r/z
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         //gc: predictionRPhi uses the cosine rule to find the phi of the 3rd point at radius, assuming the pairCurvature range [-c,+c]
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();                                    // Now recover hits in bounding box...
0362       float prmin = phiRange.min(), prmax = phiRange.max();  //get contiguous range
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         // limit error to 90 degree
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();  // they will be modified!
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         // to bookkeep the triplets and 3rd layers in triplet EDProducer
0445         if (tripletLastLayerIndex)
0446           tripletLastLayerIndex->push_back(il);
0447       }
0448     }
0449   }
0450   // std::cout << "found triplets " << result.size() << std::endl;
0451 }