TkSeedGenerator/src/FastHelix.cc

0001 #include "RecoTracker/TkSeedGenerator/interface/FastHelix.h"
0002 #include "RecoTracker/TkSeedGenerator/interface/FastLine.h"
0003
0004 void FastHelix::compute() {
0005   if (isValid() && (std::abs(tesla0) > 1e-3) && theCircle.rho() < maxRho)
0006     helixStateAtVertex();
0007   else
0008     straightLineStateAtVertex();
0009 }
0010
0011 void FastHelix::helixStateAtVertex() {
0012   // given the above rho>0.
0013   double rho = theCircle.rho();
0014   //remember (radius rho in cm):
0015   //rho =
0016   //100. * pt *
0017   //(10./(3.*MagneticField::inTesla(GlobalPoint(0., 0., 0.)).z()));
0018
0019   // pt = 0.01 * rho * (0.3*MagneticField::inTesla(GlobalPoint(0.,0.,0.)).z());
0020   double cm2GeV = 0.01 * 0.3 * tesla0;
0021   double pt = cm2GeV * rho;
0022
0023   // verify that rho is not toooo large
0024   double dcphi = ((outerHit().x() - theCircle.x0()) * (middleHit().x() - theCircle.x0()) +
0025                   (outerHit().y() - theCircle.y0()) * (middleHit().y() - theCircle.y0())) /
0026                  (rho * rho);
0027   if (std::abs(dcphi) >= 1.f) {
0028     straightLineStateAtVertex();
0029     return;
0030   }
0031
0032   GlobalPoint pMid(middleHit());
0033   GlobalPoint v(vertex());
0034
0035   // tangent in v (or the opposite...)
0036   double px = -cm2GeV * (v.y() - theCircle.y0());
0037   double py = cm2GeV * (v.x() - theCircle.x0());
0038   // check sign with scalar product
0039   if (px * (pMid.x() - v.x()) + py * (pMid.y() - v.y()) < 0.) {
0040     px = -px;
0041     py = -py;
0042   }
0043
0044   //calculate z0, pz
0045   //(z, R*phi) linear relation in a helix
0046   //with R, phi defined as radius and angle w.r.t. centre of circle
0047   //in transverse plane
0048   //pz = pT*(dz/d(R*phi)))
0049
0050   // VI 23/01/2012
0051   double dzdrphi = outerHit().z() - middleHit().z();
0052   dzdrphi /= rho * acos(dcphi);
0053   double pz = pt * dzdrphi;
0054
0055   TrackCharge q = 1;
0056   if (theCircle.x0() * py - theCircle.y0() * px < 0)
0057     q = -q;
0058   if (tesla0 < 0.)
0059     q = -q;
0060
0061   //VI
0062   if (useBasisVertex) {
0063     atVertex = GlobalTrajectoryParameters(basisVertex, GlobalVector(px, py, pz), q, bField);
0064   } else {
0065     double z_0 = middleHit().z();
0066     // assume v is before middleHit (opposite to outer)
0067     double ds = ((v.x() - theCircle.x0()) * (middleHit().x() - theCircle.x0()) +
0068                  (v.y() - theCircle.y0()) * (middleHit().y() - theCircle.y0())) /
0069                 (rho * rho);
0070     if (std::abs(ds) < 1.f) {
0071       ds = rho * acos(ds);
0072       z_0 -= ds * dzdrphi;
0073     } else {  // line????
0074       z_0 -= std::sqrt((middleHit() - v).perp2() / (outerHit() - middleHit()).perp2()) *
0075              (outerHit().z() - middleHit().z());
0076     }
0077
0078     //double z_old = -flfit.c()/flfit.n2();
0079     // std::cout << "v:xyz, z,old,new " << v << "   " << z_old << " " << z_0 << std::endl;
0080
0081     atVertex = GlobalTrajectoryParameters(GlobalPoint(v.x(), v.y(), z_0), GlobalVector(px, py, pz), q, bField);
0082   }
0083 }
0084
0085 void FastHelix::straightLineStateAtVertex() {
0086   //calculate GlobalTrajectoryParameters assuming straight line...
0087
0088   GlobalPoint pMid(middleHit());
0089   GlobalPoint v(vertex());
0090
0091   double dydx = 0.;
0092   double pt = 0., px = 0., py = 0.;
0093
0094   if (fabs(theCircle.n1()) > 0. || fabs(theCircle.n2()) > 0.)
0095     pt = maxPt;  // 10 TeV //else no pt
0096   if (fabs(theCircle.n2()) > 0.) {
0097     dydx = -theCircle.n1() / theCircle.n2();  //else px = 0
0098   }
0099   px = pt / sqrt(1. + dydx * dydx);
0100   py = px * dydx;
0101   // check sign with scalar product
0102   if (px * (pMid.x() - v.x()) + py * (pMid.y() - v.y()) < 0.) {
0103     px *= -1.;
0104     py *= -1.;
0105   }
0106
0107   //calculate z_0 and pz at vertex using weighted mean
0108   //z = z(r) = z0 + (dz/dr)*r
0109   //tan(theta) = dr/dz = (dz/dr)^-1
0110   //theta = atan(1./dzdr)
0111   //p = pt/sin(theta)
0112   //pz = p*cos(theta) = pt/tan(theta)
0113
0114   FastLine flfit(outerHit(), middleHit());
0115   double dzdr = -flfit.n1() / flfit.n2();
0116   double pz = pt * dzdr;
0117
0118   TrackCharge q = 1;
0119   //VI
0120
0121   if (useBasisVertex) {
0122     atVertex = GlobalTrajectoryParameters(basisVertex, GlobalVector(px, py, pz), q, bField);
0123   } else {
0124     double z_0 = -flfit.c() / flfit.n2();
0125     atVertex = GlobalTrajectoryParameters(GlobalPoint(v.x(), v.y(), z_0), GlobalVector(px, py, pz), q, bField);
0126   }
0127 }