RungeKutta/src/PathToPlane2Order.cc

0001 #include "PathToPlane2Order.h"
0002 #include "RKLocalFieldProvider.h"
0003 #include "FrameChanger.h"
0004 #include "TrackingTools/GeomPropagators/interface/HelixArbitraryPlaneCrossing.h"
0005
0006 #include <iostream>
0007
0008 std::pair<bool, double> PathToPlane2Order::operator()(const Plane& plane,
0009                                                       const Vector3D& pos,
0010                                                       const Vector3D& momentum,
0011                                                       double charge,
0012                                                       const PropagationDirection propDir) {
0013   // access to the field in field frame local coordinates
0014   RKLocalFieldProvider::Vector B = theField.inTesla(pos.x(), pos.y(), pos.z());
0015
0016   // Frame::GlobalVector localZ = Frame::GlobalVector( B.unit()); // local Z along field
0017   // transform field axis to global frame
0018   Frame::GlobalVector localZ = theFieldFrame->toGlobal(Frame::LocalVector(B.unit()));  // local Z along field
0019
0020   Frame::GlobalVector localY = localZ.cross(Frame::GlobalVector(1, 0, 0));
0021   if (localY.mag() < 0.1) {
0022     localY = localZ.cross(Frame::GlobalVector(0, 1, 0)).unit();
0023   } else {
0024     localY = localY.unit();
0025   }
0026   Frame::GlobalVector localX = localY.cross(localZ);
0027
0028   Frame::PositionType fpos(theFieldFrame->toGlobal(Frame::LocalPoint(pos)));
0029   Frame::RotationType frot(localX, localY, localZ);
0030   // frame in which the field is along Z
0031   Frame frame(fpos, frot);
0032
0033   //    cout << "PathToPlane2Order frame " << frame.position() << endl << frame.rotation() << endl;
0034
0035   // transform the position and direction to that frame
0036   Frame::LocalPoint localPos = frame.toLocal(fpos);  // same as LocalPoint(0,0,0)
0037
0038   //transform momentum from field frame to new frame via global frame
0039   Frame::GlobalVector gmom(theFieldFrame->toGlobal(Frame::LocalVector(momentum)));
0040   Frame::LocalVector localMom = frame.toLocal(gmom);
0041
0042   // transform the plane to the same frame
0043   Plane localPlane = FrameChanger::transformPlane(plane, frame);
0044   /*
0045      cout << "PathToPlane2Order input plane       " << plane.position() << endl
0046      << plane.rotation() << endl;
0047      cout << "PathToPlane2Order transformed plane " << localPlane->position() << endl
0048      << localPlane->rotation() << endl;
0049 */
0050   double k = 2.99792458e-3;
0051   double transverseMomentum = localMom.perp();  // transverse to the field
0052   if (!(transverseMomentum != 0)) {             // if (!(x!=0)) will trap both 0 and NaN
0053     //LogDebug("PathToPlane2Order_ZeroMomentum") << "Momentum transverse to the field is zero or Nan (" << transverseMomentum << ")\n";
0054     return std::pair<bool, double>(false, 0);
0055   }
0056   double curvature = -k * charge * B.mag() / transverseMomentum;
0057   /*
0058      cout << "PathToPlane2Order curvature " << curvature << endl;
0059      cout << "transverseMomentum " << transverseMomentum << endl;
0060      cout << "B.mag() " << B.mag() << endl;
0061      cout << "localZ " << localZ << endl;
0062      cout << "pos      " << pos << endl;
0063      cout << "momentum " << momentum << endl;
0064      cout << "localPos " << localPos << endl;
0065      cout << "localMom " << localMom << endl;
0066 */
0067   /*
0068     cout << "PathToPlane2Order: local pos " << localPos << " mom " << localMom
0069      << " curvature " << curvature << endl;
0070     cout << "PathToPlane2Order: local plane pos " << localPlane->position()
0071      << " normal " << localPlane->normalVector() << endl;
0072 */
0073   HelixArbitraryPlaneCrossing crossing(localPos.basicVector(), localMom.basicVector(), curvature, propDir);
0074   std::pair<bool, double> res = crossing.pathLength(localPlane);
0075
0076   return res;
0077 }