GeomPropagators/test/HelixPropagators_t.cpp

0001 #include "MagneticField/Engine/interface/MagneticField.h"
0002 #include "DataFormats/GeometrySurface/interface/Plane.h"
0003 #include "TrackingTools/GeomPropagators/interface/HelixBarrelPlaneCrossing2OrderLocal.h"
0004 #include "TrackingTools/GeomPropagators/interface/HelixBarrelPlaneCrossingByCircle.h"
0005 #include "TrackingTools/GeomPropagators/interface/OptimalHelixPlaneCrossing.h"
0006
0007 #include <algorithm>
0008 #include <cmath>
0009 #include <tuple>
0010
0011 #include <iostream>
0012
0013 typedef Surface::GlobalPoint GlobalPoint;
0014 typedef Surface::GlobalVector GlobalVector;
0015 typedef Surface::LocalPoint LocalPoint;
0016 typedef Surface::LocalVector LocalVector;
0017
0018 std::pair<LocalPoint, LocalVector> secondOrderAccurate(const GlobalPoint& startingPos,
0019                                                        const GlobalVector& startingDir,
0020                                                        double rho,
0021                                                        const Plane& plane) {
0022   typedef Basic2DVector<float> Vector2D;
0023
0024   // translate problem to local frame of the plane
0025   LocalPoint lPos = plane.toLocal(startingPos);
0026   LocalVector lDir = plane.toLocal(startingDir);
0027
0028   LocalVector yPrime = plane.toLocal(GlobalVector(0, 0, 1.f));
0029   float sinPhi = 0, cosPhi = 0;
0030   Vector2D pos;
0031   Vector2D dir;
0032
0033   sinPhi = yPrime.y();
0034   cosPhi = yPrime.x();
0035   pos = Vector2D(lPos.x() * cosPhi + lPos.y() * sinPhi, -lPos.x() * sinPhi + lPos.y() * cosPhi);
0036   dir = Vector2D(lDir.x() * cosPhi + lDir.y() * sinPhi, -lDir.x() * sinPhi + lDir.y() * cosPhi);
0037
0038   double d = -lPos.z();
0039   double x = pos.x() + dir.x() / lDir.z() * d - 0.5 * rho * d * d;
0040   double y = pos.y() + dir.y() / lDir.z() * d;
0041
0042   LocalPoint thePos(x * cosPhi - y * sinPhi, x * sinPhi + y * cosPhi, 0);
0043   float px = dir.x() + rho * d;
0044   LocalVector theDir(px * cosPhi - dir.y() * sinPhi, px * sinPhi + dir.y() * cosPhi, lDir.z());
0045
0046   return std::pair<LocalPoint, LocalVector>(thePos, theDir);
0047 }
0048
0049 void crossing1() {
0050   GlobalPoint startingPos(-7.79082, -47.6418, 9.18163);
0051   GlobalVector startingDir(-0.553982, -5.09198, 1.02212);
0052
0053   GlobalPoint pos(-2.95456, -48.2127, 3.1033);
0054
0055   double rho = 0.00223254;
0056
0057   Surface::RotationType rot(
0058       0.995292, 0.0969201, 0.000255868, 8.57131e-06, 0.00255196, -0.999997, -0.0969205, 0.995289, 0.00253912);
0059
0060   std::cout << rot << std::endl;
0061
0062   Plane plane(pos, rot);
0063   GlobalVector u = plane.normalVector();
0064   std::cout << "norm " << u << std::endl;
0065
0066   HelixBarrelPlaneCrossingByCircle precise(startingPos, startingDir, rho, alongMomentum);
0067   bool cross;
0068   double s;
0069   std::tie(cross, s) = precise.pathLength(plane);
0070
0071   HelixBarrelPlaneCrossing2OrderLocal crossing(startingPos, startingDir, rho, plane);
0072
0073   std::cout << plane.toLocal(GlobalPoint(precise.position(s))) << " "
0074             << plane.toLocal(GlobalVector(precise.direction(s))) << std::endl;
0075   std::cout << HelixBarrelPlaneCrossing2OrderLocal::positionOnly(startingPos, startingDir, rho, plane) << ' ';
0076   std::cout << crossing.position() << ' ' << crossing.direction() << std::endl;
0077
0078   LocalPoint thePos;
0079   LocalVector theDir;
0080   std::tie(thePos, theDir) = secondOrderAccurate(startingPos, startingDir, rho, plane);
0081
0082   std::cout << thePos << ' ' << theDir << std::endl;
0083 }
0084
0085 void crossing2() {
0086   GlobalPoint startingPos(-8.12604, -50.829, 9.82116);
0087   GlobalVector startingDir(-0.517536, -5.09581, 1.02212);
0088
0089   GlobalPoint pos(-2.96723, -51.4573, 14.8322);
0090
0091   Surface::RotationType rot(
0092       0.995041, 0.0994701, 0.000124443, 0.000108324, -0.00233467, 0.999997, 0.0994701, -0.995038, -0.00233387);
0093   std::cout << rot << std::endl;
0094
0095   Plane plane(pos, rot);
0096
0097   GlobalVector u = plane.normalVector();
0098   std::cout << "norm " << u << std::endl;
0099
0100   double rho = 0.00223254;
0101
0102   HelixBarrelPlaneCrossingByCircle precise(startingPos, startingDir, rho, alongMomentum);
0103   bool cross;
0104   double s;
0105   std::tie(cross, s) = precise.pathLength(plane);
0106
0107   std::cout << s << ' ' << precise.position(s) << ' ' << precise.direction(s) << std::endl;
0108   std::cout << plane.toLocal(GlobalPoint(precise.position(s))) << " "
0109             << plane.toLocal(GlobalVector(precise.direction(s))) << std::endl;
0110
0111   {
0112     OptimalHelixPlaneCrossing ocrossing(plane,
0113                                         HelixPlaneCrossing::PositionType(startingPos),
0114                                         HelixPlaneCrossing::DirectionType(startingDir),
0115                                         rho,
0116                                         alongMomentum);
0117     auto& crossing = *ocrossing;
0118     bool cross;
0119     double s;
0120     std::tie(cross, s) = crossing.pathLength(plane);
0121     if (!cross)
0122       std::cout << "missed!" << std::endl;
0123     else {
0124       std::cout << s << ' ' << crossing.position(s) << ' ' << crossing.direction(s) << std::endl;
0125       std::cout << plane.toLocal(GlobalPoint(crossing.position(s))) << " "
0126                 << plane.toLocal(GlobalVector(crossing.direction(s))) << std::endl;
0127     }
0128   }
0129
0130   HelixBarrelPlaneCrossing2OrderLocal crossing(startingPos, startingDir, rho, plane);
0131   std::cout << HelixBarrelPlaneCrossing2OrderLocal::positionOnly(startingPos, startingDir, rho, plane) << ' ';
0132   std::cout << crossing.position() << ' ' << crossing.direction() << std::endl;
0133
0134   LocalPoint thePos;
0135   LocalVector theDir;
0136   std::tie(thePos, theDir) = secondOrderAccurate(startingPos, startingDir, rho, plane);
0137
0138   std::cout << thePos << ' ' << theDir << std::endl;
0139 }
0140
0141 void crossing3() {
0142   GlobalPoint startingPos(-8.12604, -50.829, 9.82116);
0143   GlobalVector startingDir(-0.517536, -5.09581, 1.02212);
0144
0145   double rho = 0.00223254;
0146
0147   GlobalPoint pos(26.2991, 36.199, 100.263);
0148   Surface::RotationType rot(
0149       -0.808996, 0.587813, 1.16466e-16, 0.587813, 0.808996, -3.78444e-17, -1.16466e-16, 3.78444e-17, -1);
0150
0151   std::cout << rot << std::endl;
0152
0153   Plane plane(pos, rot);
0154   GlobalVector u = plane.normalVector();
0155   std::cout << "norm " << u << std::endl;
0156
0157   OptimalHelixPlaneCrossing ocrossing(plane,
0158                                       HelixPlaneCrossing::PositionType(startingPos),
0159                                       HelixPlaneCrossing::DirectionType(startingDir),
0160                                       rho,
0161                                       alongMomentum);
0162   auto& crossing = *ocrossing;
0163   bool cross;
0164   double s;
0165   std::tie(cross, s) = crossing.pathLength(plane);
0166   if (!cross)
0167     std::cout << "missed!" << std::endl;
0168   else {
0169     std::cout << crossing.position(s) << ' ' << crossing.direction(s) << std::endl;
0170     std::cout << plane.toLocal(GlobalPoint(crossing.position(s))) << " "
0171               << plane.toLocal(GlobalVector(crossing.direction(s))) << std::endl;
0172   }
0173
0174   LocalPoint thePos;
0175   LocalVector theDir;
0176   std::tie(thePos, theDir) = secondOrderAccurate(startingPos, startingDir, rho, plane);
0177
0178   std::cout << thePos << ' ' << theDir << std::endl;
0179 }
0180
0181 void testHelixBarrelPlaneCrossing2OrderLocal() {
0182   crossing1();
0183   std::cout << std::endl;
0184
0185   crossing2();
0186   std::cout << std::endl;
0187
0188   crossing3();
0189 }
0190
0191 int main() {
0192   testHelixBarrelPlaneCrossing2OrderLocal();
0193
0194   return 0;
0195 }