PathToPlane2Order.cc

Go to the documentation of this file.

#include "PathToPlane2Order.h"
#include "RKLocalFieldProvider.h"
#include "FrameChanger.h"
#include "TrackingTools/GeomPropagators/interface/HelixArbitraryPlaneCrossing.h"
 
#include <iostream>
 
std::pair<bool, double> PathToPlane2Order::operator()(const Plane& plane,
                                                      const Vector3D& pos,
                                                      const Vector3D& momentum,
                                                      double charge,
                                                      const PropagationDirection propDir) {
  // access to the field in field frame local coordinates
  RKLocalFieldProvider::Vector B = theField.inTesla(pos.x(), pos.y(), pos.z());
 
  // Frame::GlobalVector localZ = Frame::GlobalVector( B.unit()); // local Z along field
  // transform field axis to global frame
  Frame::GlobalVector localZ = theFieldFrame->toGlobal(Frame::LocalVector(B.unit()));  // local Z along field
 
  Frame::GlobalVector localY = localZ.cross(Frame::GlobalVector(1, 0, 0));
  if (localY.mag() < 0.1) {
    localY = localZ.cross(Frame::GlobalVector(0, 1, 0)).unit();
  } else {
    localY = localY.unit();
  }
  Frame::GlobalVector localX = localY.cross(localZ);
 
  Frame::PositionType fpos(theFieldFrame->toGlobal(Frame::LocalPoint(pos)));
  Frame::RotationType frot(localX, localY, localZ);
  // frame in which the field is along Z
  Frame frame(fpos, frot);
 
  //    cout << "PathToPlane2Order frame " << frame.position() << endl << frame.rotation() << endl;
 
  // transform the position and direction to that frame
  Frame::LocalPoint localPos = frame.toLocal(fpos);  // same as LocalPoint(0,0,0)
 
  //transform momentum from field frame to new frame via global frame
  Frame::GlobalVector gmom(theFieldFrame->toGlobal(Frame::LocalVector(momentum)));
  Frame::LocalVector localMom = frame.toLocal(gmom);
 
  // transform the plane to the same frame
  Plane localPlane = FrameChanger::transformPlane(plane, frame);
  /*
     cout << "PathToPlane2Order input plane       " << plane.position() << endl 
     << plane.rotation() << endl;
     cout << "PathToPlane2Order transformed plane " << localPlane->position() << endl 
     << localPlane->rotation() << endl;
*/
  double k = 2.99792458e-3;
  double transverseMomentum = localMom.perp();  // transverse to the field
  if (!(transverseMomentum != 0)) {             // if (!(x!=0)) will trap both 0 and NaN
    //LogDebug("PathToPlane2Order_ZeroMomentum") << "Momentum transverse to the field is zero or Nan (" << transverseMomentum << ")\n";
    return std::pair<bool, double>(false, 0);
  }
  double curvature = -k * charge * B.mag() / transverseMomentum;
  /*
     cout << "PathToPlane2Order curvature " << curvature << endl;
     cout << "transverseMomentum " << transverseMomentum << endl;
     cout << "B.mag() " << B.mag() << endl;
     cout << "localZ " << localZ << endl;
     cout << "pos      " << pos << endl;
     cout << "momentum " << momentum << endl;
     cout << "localPos " << localPos << endl;
     cout << "localMom " << localMom << endl;
*/
  /*
    cout << "PathToPlane2Order: local pos " << localPos << " mom " << localMom 
     << " curvature " << curvature << endl;
    cout << "PathToPlane2Order: local plane pos " << localPlane->position() 
     << " normal " << localPlane->normalVector() << endl;
*/
  HelixArbitraryPlaneCrossing crossing(localPos.basicVector(), localMom.basicVector(), curvature, propDir);
  std::pair<bool, double> res = crossing.pathLength(localPlane);
 
  return res;
}