ConversionSeedGenerators/interface/IdealHelixParameters.icc

0001 void IdealHelixParameters::setData(const reco::Track* track, const math::XYZPoint& ref) {
0002   setData(track, math::XYZVector(ref.x(), ref.y(), ref.z()));
0003 }
0004
0005 void IdealHelixParameters::setData(const reco::Track* track, const math::XYZVector& refPoint) {
0006   if (track != nullptr && _track != track) {
0007     _track = track;
0008     _refPoint = refPoint;
0009     calculate();
0010   }
0011
0012   if ((refPoint - _refPoint).r() != 0) {
0013     _refPoint = refPoint;
0014     calculate();
0015   }
0016 }
0017
0018 void IdealHelixParameters::calculate() {
0019   LogDebug("IdealHelixParameters") << "[IdealHelixParameters] "
0020                                    << "\n\t "
0021                                    << "refPoint \t" << _refPoint << "\n\t "
0022                                    << "innerDetid    \t" << _track->innerDetId() << "\n\t "
0023                                    << "innerMomentum \t" << _track->innerMomentum() << "\n\t "
0024                                    << "innerPosition \t" << _track->innerPosition();
0025
0026   evalCircleCenter();
0027
0028   evalTangentPoint();
0029
0030   evalMomentumatTangentPoint();
0031
0032   /* DEBUG
0033     math::XYZVector _checkMomentumAtTangentoPoint (
0034     cos_angle*innerMomentum.x()+sin_angle*innerMomentum.y(),
0035     -1*sin_angle*innerMomentum.x()+cos_angle*innerMomentum.y(),
0036     innerMomentum.z()
0037     );
0038     math::XYZVector cp_xy(_checkMomentumAtTangentoPoint.x(),_checkMomentumAtTangentoPoint.y(),0);
0039     float _checktransverseIP = (r_xy.Cross(cp_xy)).z()/cp_xy.rho();
0040
0041     ss
0042     << "\n\t "<< "_checkMomentumAtTangentoPoint \t" <<_checkMomentumAtTangentoPoint
0043     << "\n\t "<< "_checktransverseIP \t" <<_checktransverseIP;
0044   */
0045 }
0046
0047 void IdealHelixParameters::evalCircleCenter() {
0048   GlobalPoint innerPos(_track->innerPosition().x(), _track->innerPosition().y(), _track->innerPosition().z());
0049   GlobalVector innerMom(_track->innerMomentum().x(), _track->innerMomentum().y(), _track->innerMomentum().z());
0050
0051   GlobalTrajectoryParameters globalTrajParam(innerPos, innerMom, _track->charge(), _magnField);
0052
0053   _radius = 1. / fabs(globalTrajParam.transverseCurvature());
0054   float phi = innerMom.phi();
0055
0056   //eval Center of the Circumference passing in track->innserPosition
0057   float Xc = innerPos.x() + _track->charge() * _radius * sin(phi);
0058   float Yc = innerPos.y() - _track->charge() * _radius * cos(phi);
0059
0060   _circleCenter.SetXYZ(Xc, Yc, innerPos.z());  //NB Z component is useless
0061
0062   LogDebug("IdealHelixParameters") << "\n\t "
0063                                    << "circle center \t" << _circleCenter << "\n\t "
0064                                    << "radius " << _radius;
0065 }
0066
0067 void IdealHelixParameters::evalMomentumatTangentPoint() {
0068   if (_tangentPoint.r() == 0) {
0069     _MomentumAtTangentPoint = math::XYZVector(0., 0., 0.);
0070     return;
0071   }
0072
0073   math::XYZVector innerPosition(_track->innerPosition().X(), _track->innerPosition().Y(), _track->innerPosition().Z());
0074   math::XYZVector innerMomentum(_track->innerMomentum());
0075
0076   math::XYZVector pCi = innerPosition - _circleCenter;
0077   math::XYZVector pCT = _tangentPoint - _circleCenter;
0078
0079   math::XYZVector pCi_xy(pCi.x(), pCi.y(), 0);
0080   math::XYZVector pCT_xy(pCT.x(), pCT.y(), 0);
0081
0082   float cos_angle = pCi_xy.Dot(pCT_xy) / pCi_xy.rho() / pCT_xy.rho();
0083   float sin_angle = pCi_xy.Cross(pCT_xy).z() / pCi_xy.rho() / pCT_xy.rho();
0084
0085   _MomentumAtTangentPoint = math::XYZVector(cos_angle * innerMomentum.x() - sin_angle * innerMomentum.y(),
0086                                             sin_angle * innerMomentum.x() + cos_angle * innerMomentum.y(),
0087                                             innerMomentum.z());
0088
0089   math::XYZVector r_(_tangentPoint.x(), _tangentPoint.y(), 0);
0090   math::XYZVector p_(_MomentumAtTangentPoint.x(), _MomentumAtTangentPoint.y(), 0);
0091
0092   _transverseIP = (r_.x() * p_.y() - r_.y() * p_.x()) / p_.rho() / r_.rho();
0093   _rotationAngle = atan(sin_angle / cos_angle);
0094
0095   LogDebug("IdealHelixParameters") << "\n\t "
0096                                    << "_MomentumAtTangentPoint \t" << _MomentumAtTangentPoint << "\n\t "
0097                                    << "sin_angle \t" << sin_angle << " angle \t" << asin(sin_angle) << "\n\t "
0098                                    << "cos_angle \t" << cos_angle << " angle \t" << acos(cos_angle) << "\n\t "
0099                                    << "_rotationAngle \t" << _rotationAngle << "\n\t "
0100                                    << "_transverseIP \t" << _transverseIP << "\n\t "
0101                                    << " check similitude pz/pt "
0102                                    << _track->innerMomentum().z() / _track->innerMomentum().rho() << " pZ/pRho "
0103                                    << innerPosition.z() / innerPosition.rho();
0104 }
0105
0106 void IdealHelixParameters::evalTangentPoint() {
0107   math::XYZVector innerPosition(_track->innerPosition().X(), _track->innerPosition().Y(), _track->innerPosition().Z());
0108
0109   math::XYZVector vL = _circleCenter - _refPoint;
0110   float l = vL.rho();
0111
0112   if (l < _radius) {
0113     //refpoint is inside the circle
0114     _tangentPoint.SetXYZ(0., 0., 0.);
0115     return;
0116   }
0117
0118   float sin_alpha = _radius / l;
0119
0120   //there are two possible tangents, with Point of tangence T1 and T2, and same distance T from _refPoint
0121   float T = sqrt(l * l - _radius * _radius);
0122
0123   math::XYZVector vLperp(-vL.y(), vL.x(), 0);  //NB: z component not correct
0124   math::XYZVector tmpT1 = (1 - sin_alpha * sin_alpha) * vL + T / l * sin_alpha * vLperp;
0125   math::XYZVector tmpT2 = (1 - sin_alpha * sin_alpha) * vL - T / l * sin_alpha * vLperp;
0126
0127   if ((tmpT1 - innerPosition).rho() < (tmpT2 - innerPosition).rho())
0128     _tangentPoint = tmpT1 + _refPoint;
0129   else
0130     _tangentPoint = tmpT2 + _refPoint;
0131
0132   //Fix the Z component
0133   _tangentPoint.SetZ((_tangentPoint.rho() - _refPoint.rho()) *
0134                          (_track->innerMomentum().z() / _track->innerMomentum().rho()) +
0135                      _refPoint.z());
0136
0137   LogDebug("IdealHelixParameters") << "\n\t "
0138                                    << "tangent Point \t" << _tangentPoint;
0139 }
0140
0141 bool IdealHelixParameters::isTangentPointDistanceLessThan(float rMax,
0142                                                           const reco::Track* track,
0143                                                           const math::XYZVector& refPoint) {
0144   setData(track, refPoint);
0145   if (GetTangentPoint().r() == 0) {
0146     //this case means a null results on the IdealHelixParameters side
0147     return true;
0148   }
0149
0150   if (GetTangentPoint().rho() < rMax) {
0151     //this case means a track that has the tangent point nearby the primary vertex
0152     // if the track is primary, this number tends to be the primary vertex itself
0153     //Rejecting all the potential photon conversions having a "vertex" inside the beampipe
0154     //We should not miss too much, seen that the conversions at the beam pipe are the better reconstructed
0155     return true;
0156   }
0157
0158   return false;
0159 }