PixelSeeding/test/PixelTriplets_InvPrbl_t.cpp

0001 #include "RecoTracker/PixelSeeding/plugins/ThirdHitPredictionFromInvParabola.cc"
0002
0003 #include <cmath>
0004 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0005 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0006
0007 #include <iostream>
0008
0009 // oldcode
0010 template <class T>
0011 class MappedPoint {
0012 public:
0013   MappedPoint() : theU(0), theV(0), pRot(0) {}
0014   MappedPoint(const T& aU, const T& aV, const TkRotation<T>* aRot) : theU(aU), theV(aV), pRot(aRot) {}
0015   MappedPoint(const Basic2DVector<T>& point, const TkRotation<T>* aRot) : pRot(aRot) {
0016     T invRadius2 = T(1) / point.mag2();
0017     Basic3DVector<T> rotated = (*pRot) * point;
0018     theU = rotated.x() * invRadius2;
0019     theV = rotated.y() * invRadius2;
0020   }
0021   T u() const { return theU; }
0022   T v() const { return theV; }
0023   Basic2DVector<T> unmap() const {
0024     T radius2 = T(1) / (theU * theU + theV * theV);
0025     Basic3DVector<T> tmp = (*pRot).multiplyInverse(Basic2DVector<T>(theU, theV));
0026     return Basic2DVector<T>(tmp.x() * radius2, tmp.y() * radius2);
0027   }
0028
0029 private:
0030   T theU, theV;
0031   const TkRotation<T>* pRot;
0032 };
0033
0034 typedef MappedPoint<double> PointUV;
0035
0036 void oldCode(const GlobalPoint& P1, const GlobalPoint& P2) {
0037   typedef TkRotation<double> Rotation;
0038   typedef Basic2DVector<double> Point2D;
0039
0040   GlobalVector aX = GlobalVector(P1.x(), P1.y(), 0.).unit();
0041   GlobalVector aY(-aX.y(), aX.x(), 0.);
0042   GlobalVector aZ(0., 0., 1.);
0043   TkRotation<double> theRotation = Rotation(aX, aY, aZ);
0044
0045   PointUV p1(Point2D(P1.x(), P1.y()), &theRotation);
0046   PointUV p2(Point2D(P2.x(), P2.y()), &theRotation);
0047
0048   std::cout << "\nold for " << P1 << ", " << P2 << std::endl;
0049   std::cout << theRotation << std::endl;
0050   std::cout << p1.u() << " " << p1.v() << std::endl;
0051   std::cout << p2.u() << " " << p2.v() << std::endl;
0052   std::cout << p1.unmap() << std::endl;
0053   std::cout << p2.unmap() << std::endl;
0054 }
0055
0056 inline Basic2DVector<double> transform(Basic2DVector<double> const& p, TkRotation2D<double> const& theRotation) {
0057   return theRotation.rotate(p) / p.mag2();
0058 }
0059
0060 inline Basic2DVector<double> transformBack(Basic2DVector<double> const& p, TkRotation2D<double> const& theRotation) {
0061   return theRotation.rotateBack(p) / p.mag2();
0062 }
0063
0064 void newCode(const GlobalPoint& P1, const GlobalPoint& P2) {
0065   typedef TkRotation2D<double> Rotation;
0066   typedef Basic2DVector<double> Point2D;
0067
0068   Rotation theRotation = Rotation(Basic2DVector<double>(P1.basicVector().xy()));
0069   Point2D p1 = transform(Basic2DVector<double>(P1.basicVector().xy()), theRotation);  // (1./P1.xy().mag(),0);
0070   Point2D p2 = transform(Basic2DVector<double>(P2.basicVector().xy()), theRotation);
0071
0072   std::cout << "\nnew for " << P1 << ", " << P2 << std::endl;
0073   std::cout << theRotation << std::endl;
0074   std::cout << p1.x() << " " << p1.y() << std::endl;
0075   std::cout << p2.x() << " " << p2.y() << std::endl;
0076   std::cout << transformBack(p1, theRotation) << std::endl;
0077   std::cout << transformBack(p2, theRotation) << std::endl;
0078 }
0079
0080 namespace test {
0081   namespace PixelTriplets_InvPrbl_t {
0082     int test() {
0083       GlobalPoint P1(3., 4., 7.);
0084       GlobalPoint P2(-2., 5., 7.);
0085
0086       oldCode(P1, P2);
0087       newCode(P1, P2);
0088
0089       oldCode(P2, P1);
0090       newCode(P2, P1);
0091
0092       {
0093         ThirdHitPredictionFromInvParabola pred(P1, P2, 0.2, 0.05, 0.1);
0094         std::cout << "ip min, max " << pred.theIpRangePlus.min() << " " << pred.theIpRangePlus.max() << "  "
0095                   << pred.theIpRangeMinus.min() << " " << pred.theIpRangeMinus.max() << std::endl;
0096         std::cout << "A,B +pos " << pred.coeffA(0.1) << " " << pred.coeffB(0.1) << std::endl;
0097         std::cout << "A,B -pos " << pred.coeffA(-0.1) << " " << pred.coeffB(-0.1) << std::endl;
0098
0099         auto rp = pred.rangeRPhi(5., 1);
0100         auto rn = pred.rangeRPhi(5., -1);
0101         std::cout << "range " << rp.min() << " " << rp.max() << " " << rn.min() << " " << rn.max() << std::endl;
0102       }
0103
0104       ThirdHitPredictionFromInvParabola pred(-1.092805, 4.187564, -2.361283, 7.892722, 0.111413, 0.019043, 0.032000);
0105       std::cout << "ip min, max " << pred.theIpRangePlus.min() << " " << pred.theIpRangePlus.max() << "  "
0106                 << pred.theIpRangeMinus.min() << " " << pred.theIpRangeMinus.max() << std::endl;
0107       {
0108         auto rp = pred.rangeRPhi(11.4356, 1);
0109         auto rn = pred.rangeRPhi(11.4356, -1);
0110         std::cout << "range " << rp.min() << " " << rp.max() << " " << rn.min() << " " << rn.max() << std::endl;
0111       }
0112       {
0113         auto rp = pred.rangeRPhi(13.2131, 1);
0114         auto rn = pred.rangeRPhi(13.2131, -1);
0115         std::cout << "range " << rp.min() << " " << rp.max() << " " << rn.min() << " " << rn.max() << std::endl;
0116       }
0117
0118       return 0;
0119     }
0120   }  // namespace PixelTriplets_InvPrbl_t
0121 }  // namespace test
0122
0123 int main() { return test::PixelTriplets_InvPrbl_t::test(); }