TkTrackingRegions/src/RectangularEtaPhiTrackingRegion.cc

0001 #include <cmath>
0002 #include "RecoTracker/TkTrackingRegions/interface/RectangularEtaPhiTrackingRegion.h"
0003
0004 #include "OuterEstimator.h"
0005
0006 #include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h"
0007 #include "TrackingTools/TrajectoryParametrization/interface/LocalTrajectoryParameters.h"
0008 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0009 #include "TrackingTools/GeomPropagators/interface/StraightLinePropagator.h"
0010 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0011
0012 #include "TrackingTools/DetLayers/interface/BarrelDetLayer.h"
0013 #include "TrackingTools/DetLayers/interface/ForwardDetLayer.h"
0014 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0015 #include "RecoTracker/TkTrackingRegions/interface/HitRCheck.h"
0016 #include "RecoTracker/TkTrackingRegions/interface/HitZCheck.h"
0017 #include "RecoTracker/TkTrackingRegions/interface/HitEtaCheck.h"
0018 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
0019 #include "RecoTracker/TkMSParametrization/interface/PixelRecoUtilities.h"
0020 #include "RecoTracker/TkMSParametrization/interface/MultipleScatteringParametrisationMaker.h"
0021 #include "RecoTracker/TkMSParametrization/interface/MultipleScatteringParametrisation.h"
0022
0023 #include "MagneticField/Engine/interface/MagneticField.h"
0024
0025 #include "RecoTracker/MeasurementDet/interface/MeasurementTracker.h"
0026 #include "RecoTracker/MeasurementDet/interface/MeasurementTrackerEvent.h"
0027 #include "TrackingTools/MeasurementDet/interface/LayerMeasurements.h"
0028 #include "TrackingTools/PatternTools/interface/TrajectoryMeasurement.h"
0029 #include "DataFormats/GeometrySurface/interface/BoundPlane.h"
0030 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h"
0031 #include "TrackingTools/KalmanUpdators/interface/EtaPhiMeasurementEstimator.h"
0032 #include "DataFormats/Math/interface/normalizedPhi.h"
0033
0034 #include <iostream>
0035 #include <algorithm>
0036 #include <cctype>
0037
0038 namespace {
0039   template <class T>
0040   T sqr(T t) {
0041     return t * t;
0042   }
0043 }  // namespace
0044
0045 using namespace PixelRecoUtilities;
0046 using namespace std;
0047
0048 void RectangularEtaPhiTrackingRegion::checkTracks(reco::TrackCollection const& tracks, std::vector<bool>& mask) const {
0049   const math::XYZPoint regOrigin(origin().x(), origin().y(), origin().z());
0050   auto phi0 = phiDirection() + 0.5 * (phiMargin().right() - phiMargin().left());
0051   auto dphi = 0.5 * (phiMargin().right() + phiMargin().left());
0052
0053   assert(mask.size() == tracks.size());
0054   int i = -1;
0055   for (auto const& track : tracks) {
0056     i++;
0057     if (mask[i])
0058       continue;
0059
0060     if (track.pt() < ptMin()) {
0061       continue;
0062     }
0063     if (!etaRange().inside(track.eta())) {
0064       continue;
0065     }
0066     if (std::abs(proxim(track.phi(), phi0) - phi0) > dphi) {
0067       continue;
0068     }
0069     if (std::abs(track.dxy(regOrigin)) > originRBound()) {
0070       continue;
0071     }
0072     if (std::abs(track.dz(regOrigin)) > originZBound()) {
0073       continue;
0074     }
0075     mask[i] = true;
0076   }
0077 }
0078
0079 RectangularEtaPhiTrackingRegion::UseMeasurementTracker RectangularEtaPhiTrackingRegion::stringToUseMeasurementTracker(
0080     const std::string& name) {
0081   std::string tmp = name;
0082   std::transform(tmp.begin(), tmp.end(), tmp.begin(), ::tolower);
0083   if (tmp == "never")
0084     return UseMeasurementTracker::kNever;
0085   if (tmp == "forsistrips")
0086     return UseMeasurementTracker::kForSiStrips;
0087   if (tmp == "always")
0088     return UseMeasurementTracker::kAlways;
0089   throw cms::Exception("Configuration") << "Got invalid string '" << name
0090                                         << "', valid values are 'Never', 'ForSiStrips', 'Always' (case insensitive)";
0091 }
0092
0093 void RectangularEtaPhiTrackingRegion::initEtaRange(const GlobalVector& dir, const Margin& margin) {
0094   float eta = dir.eta();
0095   theEtaRange = Range(eta - margin.left(), eta + margin.right());
0096   theLambdaRange = Range(std::sinh(theEtaRange.min()), std::sinh(theEtaRange.max()));
0097   theMeanLambda = std::sinh(theEtaRange.mean());
0098 }
0099
0100 std::unique_ptr<HitRZCompatibility> RectangularEtaPhiTrackingRegion::checkRZOld(const DetLayer* layer,
0101                                                                                 const Hit& outerHit,
0102                                                                                 const DetLayer* outerlayer) const {
0103   bool isBarrel = (layer->location() == GeomDetEnumerators::barrel);
0104   GlobalPoint ohit = outerHit->globalPosition();
0105   float outerred_r = std::sqrt(sqr(ohit.x() - origin().x()) + sqr(ohit.y() - origin().y()));
0106   PixelRecoPointRZ outer(outerred_r, ohit.z());
0107
0108   float zMinOrigin = origin().z() - originZBound();
0109   float zMaxOrigin = origin().z() + originZBound();
0110
0111   if (!thePrecise) {
0112     float vcotMin = (outer.z() > zMaxOrigin) ? (outer.z() - zMaxOrigin) / (outer.r() + originRBound())
0113                                              : (outer.z() - zMaxOrigin) / (outer.r() - originRBound());
0114     float vcotMax = (outer.z() > zMinOrigin) ? (outer.z() - zMinOrigin) / (outer.r() - originRBound())
0115                                              : (outer.z() - zMinOrigin) / (outer.r() + originRBound());
0116     float cotRight = std::max(vcotMin, theLambdaRange.min());
0117     float cotLeft = std::min(vcotMax, theLambdaRange.max());
0118     return std::make_unique<HitEtaCheck>(isBarrel, outer, cotLeft, cotRight);
0119   }
0120
0121   float outerZscatt = 0;
0122   float innerScatt = 0;
0123   //CHECK
0124   if (theUseMS) {
0125     MultipleScatteringParametrisation oSigma = theMSMaker->parametrisation(layer);
0126     float cotThetaOuter = theMeanLambda;
0127     float sinThetaOuterInv = std::sqrt(1.f + sqr(cotThetaOuter));
0128     outerZscatt = 3.f * oSigma(ptMin(), cotThetaOuter) * sinThetaOuterInv;
0129   }
0130
0131   PixelRecoLineRZ boundL(outer, theLambdaRange.max());
0132   PixelRecoLineRZ boundR(outer, theLambdaRange.min());
0133   float zMinLine = boundL.zAtR(0.) - outerZscatt;
0134   float zMaxLine = boundR.zAtR(0.) + outerZscatt;
0135   PixelRecoPointRZ vtxL(0., max(zMinLine, zMinOrigin));
0136   PixelRecoPointRZ vtxR(0., min(zMaxLine, zMaxOrigin));
0137   PixelRecoPointRZ vtxMean(0., (vtxL.z() + vtxR.z()) * 0.5f);
0138   //CHECK
0139
0140   if (theUseMS) {
0141     MultipleScatteringParametrisation iSigma = theMSMaker->parametrisation(layer);
0142
0143     innerScatt =
0144         3.f * (outerlayer ? iSigma(ptMin(), vtxMean, outer, outerlayer->seqNum()) : iSigma(ptMin(), vtxMean, outer));
0145
0146     // innerScatt = 3.f *iSigma( ptMin(), vtxMean, outer);
0147   }
0148
0149   SimpleLineRZ leftLine(vtxL, outer);
0150   SimpleLineRZ rightLine(vtxR, outer);
0151
0152   HitRZConstraint rzConstraint(leftLine, rightLine);
0153   auto cotTheta = std::abs(leftLine.cotLine() + rightLine.cotLine()) * 0.5f;
0154
0155   // std::cout << "RectangularEtaPhiTrackingRegion " << outer.r()<<','<< outer.z() << " " << innerScatt << " " << cotTheta << " " <<  hitZErr <<  std::endl;
0156
0157   if (isBarrel) {
0158     auto sinThetaInv = std::sqrt(1.f + sqr(cotTheta));
0159     auto corr = innerScatt * sinThetaInv;
0160     return std::make_unique<HitZCheck>(rzConstraint, HitZCheck::Margin(corr, corr));
0161   } else {
0162     auto cosThetaInv = std::sqrt(1.f + sqr(1.f / cotTheta));
0163     auto corr = innerScatt * cosThetaInv;
0164     return std::make_unique<HitRCheck>(rzConstraint, HitRCheck::Margin(corr, corr));
0165   }
0166 }
0167
0168 std::unique_ptr<MeasurementEstimator> RectangularEtaPhiTrackingRegion::estimator(const BarrelDetLayer* layer) const {
0169   using Algo = HitZCheck;
0170
0171   // det dimensions
0172   float halfLength = 0.5f * layer->surface().bounds().length();
0173   float halfThickness = 0.5f * layer->surface().bounds().thickness();
0174   float z0 = layer->position().z();
0175   float radius = layer->specificSurface().radius();
0176
0177   // det ranges
0178   Range detRWindow(radius - halfThickness, radius + halfThickness);
0179   Range detZWindow(z0 - halfLength, z0 + halfLength);
0180
0181   // z prediction, skip if not intersection
0182   HitZCheck zPrediction(rzConstraint());
0183   Range hitZWindow = zPrediction.range(detRWindow.min()).intersection(detZWindow);
0184   if (hitZWindow.empty())
0185     return nullptr;
0186
0187   // phi prediction
0188   OuterHitPhiPrediction phiPrediction = phiWindow(*theField);
0189
0190   //
0191   // optional corrections for tolerance (mult.scatt, error, bending)
0192   //
0193   OuterHitPhiPrediction::Range phiRange;
0194   if (thePrecise) {
0195     auto invR = 1.f / radius;
0196     auto cotTheta = (hitZWindow.mean() - origin().z()) * invR;
0197     auto sinThetaInv = std::sqrt(1.f + sqr(cotTheta));
0198     MultipleScatteringParametrisation msSigma = theMSMaker->parametrisation(layer);
0199     auto scatt = 3.f * msSigma(ptMin(), cotTheta);
0200     auto bendR = longitudinalBendingCorrection(radius, ptMin(), *theField);
0201
0202     float hitErrRPhi = 0.;
0203     float hitErrZ = 0.;
0204     float corrPhi = (scatt + hitErrRPhi) * invR;
0205     float corrZ = scatt * sinThetaInv + bendR * std::abs(cotTheta) + hitErrZ;
0206
0207     phiPrediction.setTolerance(corrPhi);
0208     zPrediction.setTolerance(HitZCheck::Margin(corrZ, corrZ));
0209
0210     //
0211     // hit ranges in det
0212     //
0213     OuterHitPhiPrediction::Range phi1 = phiPrediction(detRWindow.min());
0214     OuterHitPhiPrediction::Range phi2 = phiPrediction(detRWindow.max());
0215     phiRange = Range(std::min(phi1.min(), phi2.min()), std::max(phi1.max(), phi2.max()));
0216     Range w1 = zPrediction.range(detRWindow.min());
0217     Range w2 = zPrediction.range(detRWindow.max());
0218     hitZWindow = Range(std::min(w1.min(), w2.min()), std::max(w1.max(), w2.max())).intersection(detZWindow);
0219   } else {
0220     phiRange = phiPrediction(detRWindow.mean());
0221   }
0222
0223   return std::make_unique<OuterEstimator<Algo>>(OuterDetCompatibility(layer, phiRange, detRWindow, hitZWindow),
0224                                                 OuterHitCompatibility<Algo>(phiPrediction, zPrediction));
0225 }
0226
0227 std::unique_ptr<MeasurementEstimator> RectangularEtaPhiTrackingRegion::estimator(const ForwardDetLayer* layer) const {
0228   using Algo = HitRCheck;
0229   // det dimensions, ranges
0230   float halfThickness = 0.5f * layer->surface().bounds().thickness();
0231   float zLayer = layer->position().z();
0232   Range detZWindow(zLayer - halfThickness, zLayer + halfThickness);
0233   Range detRWindow(layer->specificSurface().innerRadius(), layer->specificSurface().outerRadius());
0234
0235   // r prediction, skip if not intersection
0236   HitRCheck rPrediction(rzConstraint());
0237   Range hitRWindow = rPrediction.range(zLayer).intersection(detRWindow);
0238   if (hitRWindow.empty())
0239     return nullptr;
0240
0241   // phi prediction
0242   OuterHitPhiPrediction phiPrediction = phiWindow(*theField);
0243   OuterHitPhiPrediction::Range phiRange = phiPrediction(detRWindow.max());
0244
0245   //
0246   // optional corrections for tolerance (mult.scatt, error, bending)
0247   //
0248   if (thePrecise) {
0249     float cotTheta = (detZWindow.mean() - origin().z()) / hitRWindow.mean();
0250     float cosThetaInv = std::sqrt(1 + sqr(cotTheta)) / cotTheta;
0251     MultipleScatteringParametrisation msSigma = theMSMaker->parametrisation(layer);
0252     float scatt = 3.f * msSigma(ptMin(), cotTheta);
0253     float bendR = longitudinalBendingCorrection(hitRWindow.max(), ptMin(), *theField);
0254     float hitErrRPhi = 0.;
0255     float hitErrR = 0.;
0256     float corrPhi = (scatt + hitErrRPhi) / detRWindow.min();
0257     float corrR = scatt * std::abs(cosThetaInv) + bendR + hitErrR;
0258
0259     phiPrediction.setTolerance(corrPhi);
0260     rPrediction.setTolerance(HitRCheck::Margin(corrR, corrR));
0261
0262     //
0263     // hit ranges in det
0264     //
0265     Range w1, w2;
0266     if (zLayer > 0) {
0267       w1 = rPrediction.range(detZWindow.min());
0268       w2 = rPrediction.range(detZWindow.max());
0269     } else {
0270       w1 = rPrediction.range(detZWindow.max());
0271       w2 = rPrediction.range(detZWindow.min());
0272     }
0273     hitRWindow = Range(w1.min(), w2.max()).intersection(detRWindow);
0274   }
0275
0276   return std::make_unique<OuterEstimator<Algo>>(OuterDetCompatibility(layer, phiRange, hitRWindow, detZWindow),
0277                                                 OuterHitCompatibility<Algo>(phiPrediction, rPrediction));
0278 }
0279
0280 OuterHitPhiPrediction RectangularEtaPhiTrackingRegion::phiWindow(const MagneticField& field) const {
0281   auto phi0 = phiDirection();
0282   return OuterHitPhiPrediction(
0283       OuterHitPhiPrediction::Range(phi0 - thePhiMargin.left(), phi0 + thePhiMargin.right()),
0284       OuterHitPhiPrediction::Range(curvature(invPtRange().min(), field), curvature(invPtRange().max(), field)),
0285       originRBound());
0286 }
0287
0288 HitRZConstraint RectangularEtaPhiTrackingRegion::rzConstraint() const {
0289   HitRZConstraint::Point pLeft, pRight;
0290   float zMin = origin().z() - originZBound();
0291   float zMax = origin().z() + originZBound();
0292   float rMin = -originRBound();
0293   float rMax = originRBound();
0294   if (theEtaRange.max() > 0) {
0295     pRight = HitRZConstraint::Point(rMin, zMax);
0296   } else {
0297     pRight = HitRZConstraint::Point(rMax, zMax);
0298   }
0299   if (theEtaRange.min() > 0.) {
0300     pLeft = HitRZConstraint::Point(rMax, zMin);
0301   } else {
0302     pLeft = HitRZConstraint::Point(rMin, zMin);
0303   }
0304   return HitRZConstraint(pLeft, theLambdaRange.min(), pRight, theLambdaRange.max());
0305 }
0306
0307 TrackingRegion::Hits RectangularEtaPhiTrackingRegion::hits(const SeedingLayerSetsHits::SeedingLayer& layer) const {
0308   TrackingRegion::Hits result;
0309
0310   //ESTIMATOR
0311
0312   const DetLayer* detLayer = layer.detLayer();
0313
0314   bool measurementMethod = false;
0315   if (theMeasurementTrackerUsage == UseMeasurementTracker::kAlways)
0316     measurementMethod = true;
0317   else if (theMeasurementTrackerUsage == UseMeasurementTracker::kForSiStrips &&
0318            GeomDetEnumerators::isTrackerStrip(detLayer->subDetector()))
0319     measurementMethod = true;
0320
0321   if (measurementMethod) {
0322     const GlobalPoint vtx = origin();
0323     GlobalVector dir = direction();
0324
0325     std::unique_ptr<MeasurementEstimator> est;
0326     if ((GeomDetEnumerators::isTrackerPixel(detLayer->subDetector()) &&
0327          GeomDetEnumerators::isBarrel(detLayer->subDetector())) ||
0328         (!theUseEtaPhi && detLayer->location() == GeomDetEnumerators::barrel)) {
0329       const BarrelDetLayer& bl = dynamic_cast<const BarrelDetLayer&>(*detLayer);
0330       est = estimator(&bl);
0331     } else if ((GeomDetEnumerators::isTrackerPixel(detLayer->subDetector()) &&
0332                 GeomDetEnumerators::isEndcap(detLayer->subDetector())) ||
0333                (!theUseEtaPhi && detLayer->location() == GeomDetEnumerators::endcap)) {
0334       const ForwardDetLayer& fl = dynamic_cast<const ForwardDetLayer&>(*detLayer);
0335       est = estimator(&fl);
0336     }
0337
0338     EtaPhiMeasurementEstimator etaPhiEstimator((theEtaRange.second - theEtaRange.first) * 0.5f,
0339                                                (thePhiMargin.left() + thePhiMargin.right()) * 0.5f);
0340     MeasurementEstimator* findDetAndHits = &etaPhiEstimator;
0341     if (est) {
0342       LogDebug("RectangularEtaPhiTrackingRegion") << "use pixel specific estimator.";
0343       findDetAndHits = est.get();
0344     } else {
0345       LogDebug("RectangularEtaPhiTrackingRegion") << "use generic etat phi estimator.";
0346     }
0347
0348     // TSOS
0349     float phi = phiDirection();
0350     // std::cout << "dir " << direction().x()/direction().perp() <<','<< direction().y()/direction().perp() << " " << sin(phi) <<','<<cos(phi)<< std::endl;
0351     Surface::RotationType rot(sin(phi), -cos(phi), 0, 0, 0, -1, cos(phi), sin(phi), 0);
0352
0353     Plane::PlanePointer surface = Plane::build(GlobalPoint(0., 0., 0.), rot);
0354     //TrajectoryStateOnSurface tsos(lpar, *surface, theField);
0355
0356     FreeTrajectoryState fts(GlobalTrajectoryParameters(vtx, dir, 1, theField));
0357     TrajectoryStateOnSurface tsos(fts, *surface);
0358
0359     // propagator
0360     StraightLinePropagator prop(theField, alongMomentum);
0361
0362     LayerMeasurements lm(theMeasurementTracker->measurementTracker(), *theMeasurementTracker);
0363
0364     auto hits = lm.recHits(*detLayer, tsos, prop, *findDetAndHits);
0365
0366     result.reserve(hits.size());
0367     for (auto h : hits) {
0368       cache.emplace_back(h);
0369       result.emplace_back(h);
0370     }
0371
0372     LogDebug("RectangularEtaPhiTrackingRegion")
0373         << " found " << hits.size() << " minus one measurements on layer: " << detLayer->subDetector();
0374     // std::cout << "RectangularEtaPhiTrackingRegion" <<" found "<< meas.size()<<" minus one measurements on layer: "<<detLayer->subDetector() << std::endl;
0375
0376   } else {
0377     //
0378     // temporary solution (actually heavily used for Pixels....)
0379     //
0380     if (detLayer->location() == GeomDetEnumerators::barrel) {
0381       const BarrelDetLayer& bl = dynamic_cast<const BarrelDetLayer&>(*detLayer);
0382       auto est = estimator(&bl);
0383       if (!est)
0384         return result;
0385       using Algo = HitZCheck;
0386       auto const& hitComp = (reinterpret_cast<OuterEstimator<Algo> const&>(*est)).hitCompatibility();
0387       auto layerHits = layer.hits();
0388       result.reserve(layerHits.size());
0389       for (auto&& ih : layerHits) {
0390         if (hitComp(*ih))
0391           result.emplace_back(ih);
0392       }
0393
0394     } else {
0395       const ForwardDetLayer& fl = dynamic_cast<const ForwardDetLayer&>(*detLayer);
0396       auto est = estimator(&fl);
0397       if (!est)
0398         return result;
0399       using Algo = HitRCheck;
0400       auto const& hitComp = (reinterpret_cast<OuterEstimator<Algo> const&>(*est)).hitCompatibility();
0401       auto layerHits = layer.hits();
0402       result.reserve(layerHits.size());
0403       for (auto&& ih : layerHits) {
0404         if (hitComp(*ih))
0405           result.emplace_back(ih);
0406       }
0407     }
0408   }
0409
0410   // std::cout << "RectangularEtaPhiTrackingRegion hits "  << result.size() << std::endl;
0411
0412   return result;
0413 }
0414
0415 std::string RectangularEtaPhiTrackingRegion::print() const {
0416   std::ostringstream str;
0417   str << TrackingRegionBase::print() << " eta: " << theEtaRange << " phi:" << thePhiMargin << "precise: " << thePrecise;
0418   return str.str();
0419 }