GeomPropagators/src/HelixBarrelPlaneCrossing2OrderLocal.cc

0001 #include "TrackingTools/GeomPropagators/interface/HelixBarrelPlaneCrossing2OrderLocal.h"
0002 #include <algorithm>
0003 #include <cmath>
0004
0005 HelixBarrelPlaneCrossing2OrderLocal::HelixBarrelPlaneCrossing2OrderLocal(const GlobalPoint& startingPos,
0006                                                                          const GlobalVector& startingDir,
0007                                                                          float rho,
0008                                                                          const Plane& plane) {
0009   typedef Basic2DVector<float> Vector2D;
0010
0011   // translate problem to local frame of the plane
0012   Plane::ToLocal toLocal(plane);
0013   LocalPoint lPos = toLocal.toLocal(startingPos);
0014   LocalVector lDir = toLocal.toLocal(startingDir);
0015
0016   // check if local frame is already special (local Y axis == global Z axis)
0017   LocalVector yPrime = toLocal.toLocal(GlobalVector(0, 0, 1.f));
0018   // LocalVector diff = yPrime - LocalVector(0,-1.f,0);
0019   float sinPhi = 0, cosPhi = 0;
0020   bool rotated;
0021   Vector2D pos;
0022   Vector2D dir;
0023
0024   if (std::abs(yPrime.y()) > 0.9999f) {
0025     // cout << "Plane already oriented, yPrime = " << yPrime << endl;
0026
0027     // we are already in the special orientation
0028     pos = Vector2D(lPos.x(), lPos.y());
0029     dir = Vector2D(lDir.x(), lDir.y());
0030     ;
0031     rotated = false;
0032   } else {
0033     // cout << "Plane needs rotation, yPrime = " << yPrime << endl;
0034
0035     // we need to rotate the problem
0036     sinPhi = yPrime.y();
0037     cosPhi = yPrime.x();
0038     pos = Vector2D(lPos.x() * cosPhi + lPos.y() * sinPhi, -lPos.x() * sinPhi + lPos.y() * cosPhi);
0039     dir = Vector2D(lDir.x() * cosPhi + lDir.y() * sinPhi, -lDir.x() * sinPhi + lDir.y() * cosPhi);
0040     rotated = true;
0041   }
0042
0043   float d = -lPos.z();
0044   float x = pos.x() + dir.x() / lDir.z() * d - 0.5f * rho * d * d;
0045   float y = pos.y() + dir.y() / lDir.z() * d;
0046
0047   //    cout << "d= " << d << ", pos.x()= " << pos.x()
0048   //         << ", dir.x()/lDir.z()= " << dir.x()/lDir.z()
0049   //         << ", 0.5*rho*d*d= " << 0.5*rho*d*d << endl;
0050
0051   if (!rotated) {
0052     thePos = LocalPoint(x, y, 0);
0053     theDir = LocalVector(dir.x() + rho * d, dir.y(), lDir.z());
0054   } else {
0055     thePos = LocalPoint(x * cosPhi - y * sinPhi, x * sinPhi + y * cosPhi, 0);
0056     float px = dir.x() + rho * d;
0057     theDir = LocalVector(px * cosPhi - dir.y() * sinPhi, px * sinPhi + dir.y() * cosPhi, lDir.z());
0058   }
0059 }
0060
0061 LocalPoint HelixBarrelPlaneCrossing2OrderLocal::positionOnly(const GlobalPoint& startingPos,
0062                                                              const GlobalVector& startingDir,
0063                                                              float rho,
0064                                                              const Plane& plane) {
0065   typedef Basic2DVector<float> Vector2D;
0066
0067   // translate problem to local frame of the plane
0068   Plane::ToLocal toLocal(plane);
0069   LocalPoint lPos = toLocal.toLocal(startingPos);
0070   LocalVector lDir = toLocal.toLocal(startingDir);
0071
0072   // check if local frame is already special (local Y axis == global Z axis)
0073   LocalVector yPrime = toLocal.toLocal(GlobalVector(0, 0, 1.f));
0074   // LocalVector diff = yPrime - LocalVector(0,-1.f,0);
0075   float sinPhi = 0, cosPhi = 0;
0076   bool rotated;
0077   Vector2D pos;
0078   Vector2D dir;
0079
0080   if (std::abs(yPrime.y()) > 0.9999f) {
0081     // cout << "Plane already oriented, yPrime = " << yPrime << endl;
0082
0083     // we are already in the special orientation
0084     pos = Vector2D(lPos.x(), lPos.y());
0085     dir = Vector2D(lDir.x(), lDir.y());
0086     ;
0087     rotated = false;
0088   } else {
0089     // cout << "Plane needs rotation, yPrime = " << yPrime << endl;
0090
0091     // we need to rotate the problem
0092     sinPhi = yPrime.y();
0093     cosPhi = yPrime.x();
0094     pos = Vector2D(lPos.x() * cosPhi + lPos.y() * sinPhi, -lPos.x() * sinPhi + lPos.y() * cosPhi);
0095     dir = Vector2D(lDir.x() * cosPhi + lDir.y() * sinPhi, -lDir.x() * sinPhi + lDir.y() * cosPhi);
0096     rotated = true;
0097   }
0098
0099   float d = -lPos.z();
0100   float x = pos.x() + dir.x() / lDir.z() * d - 0.5f * rho * d * d;
0101   float y = pos.y() + dir.y() / lDir.z() * d;
0102
0103   //    cout << "d= " << d << ", pos.x()= " << pos.x()
0104   //         << ", dir.x()/lDir.z()= " << dir.x()/lDir.z()
0105   //         << ", 0.5*rho*d*d= " << 0.5*rho*d*d << endl;
0106
0107   return (!rotated) ? LocalPoint(x, y, 0) : LocalPoint(x * cosPhi - y * sinPhi, x * sinPhi + y * cosPhi, 0);
0108 }