File indexing completed on 2023-04-15 01:47:36
0001 #ifndef RecoTracker_PixelSeeding_plugins_ThirdHitPredictionFromInvParabola_h
0002 #define RecoTracker_PixelSeeding_plugins_ThirdHitPredictionFromInvParabola_h
0003
0004
0005
0006
0007
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
0024 namespace test {
0025 namespace PixelTriplets_InvPrbl_prec {
0026 int test();
0027 }
0028 namespace PixelTriplets_InvPrbl_t {
0029 int test();
0030 }
0031 }
0032
0033 class ThirdHitPredictionFromInvParabola {
0034
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
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;
0060
0061 Range rangeRPhi(Scalar radius) const;
0062
0063
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
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
0117
0118
0119 Scalar A = coeffA(ip);
0120 Scalar B = coeffB(ip);
0121
0122
0123 Scalar ipOverR = ip / r;
0124
0125 Scalar a = 0.5 * B + ipOverR;
0126 Scalar c = -B + A * r - ipOverR;
0127
0128 Scalar delta = 1 - 4 * a * c;
0129
0130 Scalar sqrtdelta = std::sqrt(delta);
0131
0132 Scalar alpha = (-2 * c) / (1 + sqrtdelta);
0133
0134 v = alpha;
0135 Scalar d2 = 1. - v * v;
0136
0137 u = std::sqrt(d2);
0138
0139
0140 }
0141
0142 #endif