PatternTools/src/ClosestApproachInRPhi.cc

0001 #include "TrackingTools/PatternTools/interface/ClosestApproachInRPhi.h"
0002 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0003 #include "MagneticField/Engine/interface/MagneticField.h"
0004 #include "FWCore/Utilities/interface/Exception.h"
0005
0006 using namespace std;
0007
0008 bool ClosestApproachInRPhi::calculate(const TrajectoryStateOnSurface& sta, const TrajectoryStateOnSurface& stb) {
0009   TrackCharge chargeA = sta.charge();
0010   TrackCharge chargeB = stb.charge();
0011   GlobalVector momentumA = sta.globalMomentum();
0012   GlobalVector momentumB = stb.globalMomentum();
0013   GlobalPoint positionA = sta.globalPosition();
0014   GlobalPoint positionB = stb.globalPosition();
0015   paramA = sta.globalParameters();
0016   paramB = stb.globalParameters();
0017   // compute magnetic field ONCE
0018   bz = sta.freeState()->parameters().magneticField().inTesla(positionA).z() * 2.99792458e-3;
0019
0020   return compute(chargeA, momentumA, positionA, chargeB, momentumB, positionB);
0021 }
0022
0023 bool ClosestApproachInRPhi::calculate(const FreeTrajectoryState& sta, const FreeTrajectoryState& stb) {
0024   TrackCharge chargeA = sta.charge();
0025   TrackCharge chargeB = stb.charge();
0026   GlobalVector momentumA = sta.momentum();
0027   GlobalVector momentumB = stb.momentum();
0028   GlobalPoint positionA = sta.position();
0029   GlobalPoint positionB = stb.position();
0030   paramA = sta.parameters();
0031   paramB = stb.parameters();
0032   // compute magnetic field ONCE
0033   bz = sta.parameters().magneticField().inTesla(positionA).z() * 2.99792458e-3;
0034
0035   return compute(chargeA, momentumA, positionA, chargeB, momentumB, positionB);
0036 }
0037
0038 pair<GlobalPoint, GlobalPoint> ClosestApproachInRPhi::points() const {
0039   if (!status_)
0040     throw cms::Exception(
0041         "TrackingTools/PatternTools",
0042         "ClosestApproachInRPhi::could not compute track crossing. Check status before calling this method!");
0043   return pair<GlobalPoint, GlobalPoint>(posA, posB);
0044 }
0045
0046 GlobalPoint ClosestApproachInRPhi::crossingPoint() const {
0047   if (!status_)
0048     throw cms::Exception(
0049         "TrackingTools/PatternTools",
0050         "ClosestApproachInRPhi::could not compute track crossing. Check status before calling this method!");
0051   return GlobalPoint(0.5 * (posA.basicVector() + posB.basicVector()));
0052 }
0053
0054 float ClosestApproachInRPhi::distance() const {
0055   if (!status_)
0056     throw cms::Exception(
0057         "TrackingTools/PatternTools",
0058         "ClosestApproachInRPhi::could not compute track crossing. Check status before calling this method!");
0059   return (posB - posA).mag();
0060 }
0061
0062 bool ClosestApproachInRPhi::compute(const TrackCharge& chargeA,
0063                                     const GlobalVector& momentumA,
0064                                     const GlobalPoint& positionA,
0065                                     const TrackCharge& chargeB,
0066                                     const GlobalVector& momentumB,
0067                                     const GlobalPoint& positionB) {
0068   // centres and radii of track circles
0069   double xca, yca, ra;
0070   circleParameters(chargeA, momentumA, positionA, xca, yca, ra, bz);
0071   double xcb, ycb, rb;
0072   circleParameters(chargeB, momentumB, positionB, xcb, ycb, rb, bz);
0073
0074   // points of closest approach in transverse plane
0075   double xg1, yg1, xg2, yg2;
0076   int flag = transverseCoord(xca, yca, ra, xcb, ycb, rb, xg1, yg1, xg2, yg2);
0077   if (flag == 0) {
0078     status_ = false;
0079     return false;
0080   }
0081
0082   double xga, yga, zga, xgb, ygb, zgb;
0083
0084   if (flag == 1) {
0085     // two crossing points on each track in transverse plane
0086     // select point for which z-coordinates on the 2 tracks are the closest
0087     double za1 = zCoord(momentumA, positionA, ra, xca, yca, xg1, yg1);
0088     double zb1 = zCoord(momentumB, positionB, rb, xcb, ycb, xg1, yg1);
0089     double za2 = zCoord(momentumA, positionA, ra, xca, yca, xg2, yg2);
0090     double zb2 = zCoord(momentumB, positionB, rb, xcb, ycb, xg2, yg2);
0091
0092     if (abs(zb1 - za1) < abs(zb2 - za2)) {
0093       xga = xg1;
0094       yga = yg1;
0095       zga = za1;
0096       zgb = zb1;
0097     } else {
0098       xga = xg2;
0099       yga = yg2;
0100       zga = za2;
0101       zgb = zb2;
0102     }
0103     xgb = xga;
0104     ygb = yga;
0105   } else {
0106     // one point of closest approach on each track in transverse plane
0107     xga = xg1;
0108     yga = yg1;
0109     zga = zCoord(momentumA, positionA, ra, xca, yca, xga, yga);
0110     xgb = xg2;
0111     ygb = yg2;
0112     zgb = zCoord(momentumB, positionB, rb, xcb, ycb, xgb, ygb);
0113   }
0114
0115   posA = GlobalPoint(xga, yga, zga);
0116   posB = GlobalPoint(xgb, ygb, zgb);
0117   status_ = true;
0118   return true;
0119 }
0120
0121 pair<GlobalTrajectoryParameters, GlobalTrajectoryParameters> ClosestApproachInRPhi::trajectoryParameters() const {
0122   if (!status_)
0123     throw cms::Exception(
0124         "TrackingTools/PatternTools",
0125         "ClosestApproachInRPhi::could not compute track crossing. Check status before calling this method!");
0126   pair<GlobalTrajectoryParameters, GlobalTrajectoryParameters> ret(newTrajectory(posA, paramA, bz),
0127                                                                    newTrajectory(posB, paramB, bz));
0128   return ret;
0129 }
0130
0131 GlobalTrajectoryParameters ClosestApproachInRPhi::newTrajectory(const GlobalPoint& newpt,
0132                                                                 const GlobalTrajectoryParameters& oldgtp,
0133                                                                 double bz) {
0134   // First we need the centers of the circles.
0135   double qob = oldgtp.charge() / bz;
0136   double xc = oldgtp.position().x() + qob * oldgtp.momentum().y();
0137   double yc = oldgtp.position().y() - qob * oldgtp.momentum().x();
0138
0139   // and of course....
0140   double npx = (newpt.y() - yc) * (bz / oldgtp.charge());
0141   double npy = (xc - newpt.x()) * (bz / oldgtp.charge());
0142
0143   /*
0144    * old code: slow and wrong
0145    *
0146   // now we do a translation, move the center of circle to (0,0,0).
0147   double dx1 = oldgtp.position().x() - xc;
0148   double dy1 = oldgtp.position().y() - yc;
0149   double dx2 = newpt.x() - xc;
0150   double dy2 = newpt.y() - yc;
0151
0152   // now for the angles:
0153   double cosphi = ( dx1 * dx2 + dy1 * dy2 ) /
0154     ( sqrt ( dx1 * dx1 + dy1 * dy1 ) * sqrt ( dx2 * dx2 + dy2 * dy2 ));
0155   double sinphi = - oldgtp.charge() * sqrt ( 1 - cosphi * cosphi );
0156
0157   // Finally, the new momenta:
0158   double px = cosphi * oldgtp.momentum().x() - sinphi * oldgtp.momentum().y();
0159   double py = sinphi * oldgtp.momentum().x() + cosphi * oldgtp.momentum().y();
0160
0161   std::cout << px-npx << " " << py-npy << ", " << oldgtp.charge() << std::endl;
0162   */
0163
0164   GlobalVector vta(npx, npy, oldgtp.momentum().z());
0165   GlobalTrajectoryParameters gta(newpt, vta, oldgtp.charge(), &(oldgtp.magneticField()));
0166   return gta;
0167 }
0168
0169 void ClosestApproachInRPhi::circleParameters(const TrackCharge& charge,
0170                                              const GlobalVector& momentum,
0171                                              const GlobalPoint& position,
0172                                              double& xc,
0173                                              double& yc,
0174                                              double& r,
0175                                              double bz) {
0176   // compute radius of circle
0177   /** temporary code, to be replaced by call to curvature() when bug
0178    *  is fixed.
0179    */
0180   //   double bz = MagneticField::inInverseGeV(position).z();
0181
0182   // signed_r directed towards circle center, along F_Lorentz = q*v X B
0183   double qob = charge / bz;
0184   double signed_r = qob * momentum.transverse();
0185   r = abs(signed_r);
0186   /** end of temporary code
0187    */
0188
0189   // compute centre of circle
0190   // double phi = momentum.phi();
0191   // xc = signed_r*sin(phi) + position.x();
0192   // yc = -signed_r*cos(phi) + position.y();
0193   xc = position.x() + qob * momentum.y();
0194   yc = position.y() - qob * momentum.x();
0195 }
0196
0197 int ClosestApproachInRPhi::transverseCoord(double cxa,
0198                                            double cya,
0199                                            double ra,
0200                                            double cxb,
0201                                            double cyb,
0202                                            double rb,
0203                                            double& xg1,
0204                                            double& yg1,
0205                                            double& xg2,
0206                                            double& yg2) {
0207   int flag = 0;
0208   double x1, y1, x2, y2;
0209
0210   // new reference frame with origin in (cxa, cya) and x-axis
0211   // directed from (cxa, cya) to (cxb, cyb)
0212
0213   double d_ab = sqrt((cxb - cxa) * (cxb - cxa) + (cyb - cya) * (cyb - cya));
0214   if (d_ab == 0) {  // concentric circles
0215     return 0;
0216   }
0217   // elements of rotation matrix
0218   double u = (cxb - cxa) / d_ab;
0219   double v = (cyb - cya) / d_ab;
0220
0221   // conditions for circle intersection
0222   if (d_ab <= ra + rb && d_ab >= abs(rb - ra)) {
0223     // circles cross each other
0224     flag = 1;
0225
0226     // triangle (ra, rb, d_ab)
0227     double cosphi = (ra * ra - rb * rb + d_ab * d_ab) / (2 * ra * d_ab);
0228     double sinphi2 = 1. - cosphi * cosphi;
0229     if (sinphi2 < 0.) {
0230       sinphi2 = 0.;
0231       cosphi = 1.;
0232     }
0233
0234     // intersection points in new frame
0235     double sinphi = sqrt(sinphi2);
0236     x1 = ra * cosphi;
0237     y1 = ra * sinphi;
0238     x2 = x1;
0239     y2 = -y1;
0240   } else if (d_ab > ra + rb) {
0241     // circles are external to each other
0242     flag = 2;
0243
0244     // points of closest approach in new frame
0245     // are on line between 2 centers
0246     x1 = ra;
0247     y1 = 0;
0248     x2 = d_ab - rb;
0249     y2 = 0;
0250   } else if (d_ab < abs(rb - ra)) {
0251     // circles are inside each other
0252     flag = 2;
0253
0254     // points of closest approach in new frame are on line between 2 centers
0255     // choose 2 closest points
0256     double sign = 1.;
0257     if (ra <= rb)
0258       sign = -1.;
0259     x1 = sign * ra;
0260     y1 = 0;
0261     x2 = d_ab + sign * rb;
0262     y2 = 0;
0263   } else {
0264     return 0;
0265   }
0266
0267   // intersection points in global frame, transverse plane
0268   xg1 = u * x1 - v * y1 + cxa;
0269   yg1 = v * x1 + u * y1 + cya;
0270   xg2 = u * x2 - v * y2 + cxa;
0271   yg2 = v * x2 + u * y2 + cya;
0272
0273   return flag;
0274 }
0275
0276 double ClosestApproachInRPhi::zCoord(
0277     const GlobalVector& mom, const GlobalPoint& pos, double r, double xc, double yc, double xg, double yg) {
0278   // starting point
0279   double x = pos.x();
0280   double y = pos.y();
0281   double z = pos.z();
0282
0283   double px = mom.x();
0284   double py = mom.y();
0285   double pz = mom.z();
0286
0287   // rotation angle phi from starting point to crossing point (absolute value)
0288   // -- compute sin(phi/2) if phi smaller than pi/4,
0289   // -- cos(phi) if phi larger than pi/4
0290   double phi = 0.;
0291   double sinHalfPhi = sqrt((x - xg) * (x - xg) + (y - yg) * (y - yg)) / (2 * r);
0292   if (sinHalfPhi < 0.383) {  // sin(pi/8)
0293     phi = 2 * asin(sinHalfPhi);
0294   } else {
0295     double cosPhi = ((x - xc) * (xg - xc) + (y - yc) * (yg - yc)) / (r * r);
0296     if (std::abs(cosPhi) > 1)
0297       cosPhi = (cosPhi > 0 ? 1 : -1);
0298     phi = abs(acos(cosPhi));
0299   }
0300   // -- sign of phi
0301   double signPhi = ((x - xc) * (yg - yc) - (xg - xc) * (y - yc) > 0) ? 1. : -1.;
0302
0303   // sign of track angular momentum
0304   // if rotation is along angular momentum, delta z is along pz
0305   double signOmega = ((x - xc) * py - (y - yc) * px > 0) ? 1. : -1.;
0306
0307   // delta z
0308   // -- |dz| = |cos(theta) * path along helix|
0309   //         = |cos(theta) * arc length along circle / sin(theta)|
0310   double dz = signPhi * signOmega * (pz / mom.transverse()) * phi * r;
0311
0312   return z + dz;
0313 }