EgammaPhotonAlgos/src/TangentApproachInRPhi.cc

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