PixelSeeding/plugins/ThirdHitPredictionFromInvParabola.h

0001 #ifndef RecoTracker_PixelSeeding_plugins_ThirdHitPredictionFromInvParabola_h
0002 #define RecoTracker_PixelSeeding_plugins_ThirdHitPredictionFromInvParabola_h
0003
0004 /** Check phi compatibility of a hit with a track
0005     constrained by a hit pair and TrackingRegion kinematical constraints.
0006     The "inverse parabola method" is used.
0007     M.Hansroul, H.Jeremie, D.Savard NIM A270 (1998) 490.
0008  */
0009
0010 #include "DataFormats/GeometryVector/interface/Basic2DVector.h"
0011 #include "DataFormats/GeometryVector/interface/Basic3DVector.h"
0012 #include "DataFormats/GeometrySurface/interface/TkRotation.h"
0013
0014 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0015 #include "RecoTracker/TkMSParametrization/interface/PixelRecoRange.h"
0016 #include <array>
0017
0018 #include "FWCore/Utilities/interface/Visibility.h"
0019
0020 class TrackingRegion;
0021 class OrderedHitPair;
0022
0023 // Function for testing ThirdHitPredictionFromInvParabola
0024 namespace test {
0025   namespace PixelTriplets_InvPrbl_prec {
0026     int test();
0027   }
0028   namespace PixelTriplets_InvPrbl_t {
0029     int test();
0030   }
0031 }  // namespace test
0032
0033 class ThirdHitPredictionFromInvParabola {
0034   // For tests
0035   friend int test::PixelTriplets_InvPrbl_prec::test();
0036   friend int test::PixelTriplets_InvPrbl_t::test();
0037
0038 public:
0039   using Scalar = double;
0040   typedef TkRotation2D<Scalar> Rotation;
0041   typedef PixelRecoRange<float> Range;
0042   typedef PixelRecoRange<Scalar> RangeD;
0043   typedef Basic2DVector<Scalar> Point2D;
0044
0045   ThirdHitPredictionFromInvParabola() {}
0046   ThirdHitPredictionFromInvParabola(Scalar x1, Scalar y1, Scalar x2, Scalar y2, Scalar ip, Scalar curv, Scalar tolerance)
0047       : theTolerance(tolerance) {
0048     init(x1, y1, x2, y2, ip, std::abs(curv));
0049   }
0050
0051   ThirdHitPredictionFromInvParabola(
0052       const GlobalPoint &P1, const GlobalPoint &P2, Scalar ip, Scalar curv, Scalar tolerance);
0053
0054   //  inline Range operator()(Scalar radius, int charge) const { return rangeRPhiSlow(radius,charge,1); }
0055   inline Range operator()(Scalar radius, int charge) const { return rangeRPhi(radius, charge); }
0056
0057   inline Range operator()(Scalar radius) const { return rangeRPhi(radius); }
0058
0059   Range rangeRPhi(Scalar radius, int charge) const;  //  __attribute__ ((optimize(3, "fast-math")));
0060
0061   Range rangeRPhi(Scalar radius) const;
0062
0063   // Range rangeRPhiSlow(Scalar radius, int charge, int nIter=5) const;
0064
0065   void init(const GlobalPoint &P1, const GlobalPoint &P2, Scalar ip, Scalar curv) {
0066     init(P1.x(), P1.y(), P2.x(), P2.y(), ip, curv);
0067   }
0068   void init(Scalar x1, Scalar y1, Scalar x2, Scalar y2, Scalar ip, Scalar curv);
0069
0070 private:
0071   inline Scalar coeffA(Scalar impactParameter) const;
0072   inline Scalar coeffB(Scalar impactParameter) const;
0073   inline Scalar predV(Scalar u, Scalar ip) const;
0074   inline Scalar ipFromCurvature(Scalar curvature, bool pos) const;
0075
0076   Point2D transform(Point2D const &p) const { return theRotation.rotate(p) / p.mag2(); }
0077
0078   Point2D transformBack(Point2D const &p) const { return theRotation.rotateBack(p) / p.mag2(); }
0079
0080 private:
0081   Rotation theRotation;
0082   Scalar u1u2, overDu, pv, dv, su;
0083
0084   // formula is symmetric for (ip,pos) -> (-ip,neg)
0085   inline void findPointAtCurve(Scalar radius, Scalar ip, Scalar &u, Scalar &v) const;
0086
0087   RangeD theIpRangePlus, theIpRangeMinus;
0088   Scalar theTolerance;
0089   bool emptyP, emptyM;
0090 };
0091
0092 ThirdHitPredictionFromInvParabola::Scalar ThirdHitPredictionFromInvParabola::coeffA(Scalar impactParameter) const {
0093   auto c = -pv * overDu;
0094   return c - u1u2 * impactParameter;
0095 }
0096
0097 ThirdHitPredictionFromInvParabola::Scalar ThirdHitPredictionFromInvParabola::coeffB(Scalar impactParameter) const {
0098   auto c = dv * overDu;
0099   return c - su * impactParameter;
0100 }
0101
0102 ThirdHitPredictionFromInvParabola::Scalar ThirdHitPredictionFromInvParabola::ipFromCurvature(Scalar curvature,
0103                                                                                              bool pos) const {
0104   Scalar overU1u2 = 1. / u1u2;
0105   Scalar inInf = -pv * overDu * overU1u2;
0106   return (pos ? inInf : -inInf) - curvature * overU1u2 * 0.5;
0107 }
0108
0109 ThirdHitPredictionFromInvParabola::Scalar ThirdHitPredictionFromInvParabola::predV(Scalar u, Scalar ip) const {
0110   auto c = -(coeffA(ip) - coeffB(ip * u) - ip * u * u);
0111   return c;
0112 }
0113
0114 void ThirdHitPredictionFromInvParabola::findPointAtCurve(Scalar r, Scalar ip, Scalar &u, Scalar &v) const {
0115   //
0116   // assume u=(1-alpha^2/2)/r v=alpha/r
0117   // solve qudratic equation neglecting aplha^4 term
0118   //
0119   Scalar A = coeffA(ip);
0120   Scalar B = coeffB(ip);
0121
0122   // Scalar overR = 1./r;
0123   Scalar ipOverR = ip / r;  // *overR;
0124
0125   Scalar a = 0.5 * B + ipOverR;
0126   Scalar c = -B + A * r - ipOverR;
0127
0128   Scalar delta = 1 - 4 * a * c;
0129   // Scalar sqrtdelta = (delta > 0) ? std::sqrt(delta) : 0.;
0130   Scalar sqrtdelta = std::sqrt(delta);
0131   //  Scalar alpha = (-1+sqrtdelta)/(2*a);
0132   Scalar alpha = (-2 * c) / (1 + sqrtdelta);
0133
0134   v = alpha;               // *overR
0135   Scalar d2 = 1. - v * v;  // overR*overR - v*v
0136   // u = (d2 > 0) ? std::sqrt(d2) : 0.;
0137   u = std::sqrt(d2);
0138
0139   // u,v not rotated! not multiplied by 1/r
0140 }
0141
0142 #endif