PixelSeeding/plugins/PixelTripletHLTGenerator.cc

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();  // until we have moved SeedComparitor too to EDProducers
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();  //try to take account of non-centrality (?)
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;       // re-used throughout
0132   std::vector<unsigned int> foundNodes;  // re-used thoughout
0133   foundNodes.reserve(100);
0134
0135   declareDynArray(KDTreeLinkerAlgo<unsigned int>, nThirdLayers, hitTree);
0136   std::vector<float> rzError(nThirdLayers);  //save maximum errors
0137
0138   const float maxDelphi = region.ptMin() < 0.3f ? float(M_PI) / 4.f : float(M_PI) / 8.f;  // FIXME move to config??
0139   const float maxphi = M_PI + maxDelphi, minphi = -maxphi;  // increase to cater for any range
0140   const float safePhi = M_PI - maxDelphi;                   // sideband
0141
0142   const auto& field = es.getData(fieldToken_);
0143   const MultipleScatteringParametrisationMaker* msmaker = nullptr;
0144   if (useMScat) {
0145     msmaker = &es.getData(msmakerToken_);
0146   }
0147
0148   // fill the prediction vector
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;  // Initialise to extreme values in case no hits
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       //use (phi,r) for endcaps rather than (phi,z)
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);  // save it
0176       // populate side-bands
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);  // declare our bounds
0183     //add fudge factors in case only one hit and also for floating-point inaccuracy
0184     hitTree[il].build(layerTree, phiZ);  // make KDtree
0185     rzError[il] = maxErr;                //save error
0186     // std::cout << "layer " << thirdLayerDetLayer[il]->seqNum() << " " << layerTree.size() << std::endl;
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     // printf("++Constr %f %f %f %f %f %f %f\n",xi,yi,xo,yo,imppartmp,curv,extraHitRPhitolerance);
0219
0220     // std::cout << ip << ": " << point1.r() << ","<< point1.z() << " "
0221     //                        << point2.r() << ","<< point2.z() <<std::endl;
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;  // Don't bother if no hits
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         // std::cout << "++R " << radius.min() << " " << radius.max()  << std::endl;
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;  // std::sqrt(12.f); // ...and continue as before
0295       constexpr float nSigmaPhi = 3.f;
0296
0297       foundNodes.clear();  // Now recover hits in bounding box...
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       // std::cout << ip << ": " << thirdLayerDetLayer[il]->seqNum() << " " << foundNodes.size() << " " << prmin << " " << prmax << std::endl;
0326
0327       // int kk=0;
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         //if ((kk++)%100==0)
0342         //std::cout << kk << ": " << p3_u << " " << p3_v << " " << p3_phi << std::endl;
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         // limit error to 90 degree
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;  // range is empty
0362           correction.correctRPhiRange(rangeRPhi);
0363           if (checkPhiInRange(p3_phi, rangeRPhi.first * ir - phiErr, rangeRPhi.second * ir + phiErr, maxPhiErr)) {
0364             // insert here check with comparitor
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               // to bookkeep the triplets and 3rd layers in triplet EDProducer
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   // std::cout << "triplets " << result.size() << std::endl;
0385 }