File indexing completed on 2024-04-06 12:31:35
0001 #include "TrackingTools/PatternTools/interface/TwoTrackMinimumDistanceHelixHelix.h"
0002 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0003 #include "MagneticField/Engine/interface/MagneticField.h"
0004 #include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 #include "FWCore/Utilities/interface/isFinite.h"
0007
0008 using namespace std;
0009
0010 TwoTrackMinimumDistanceHelixHelix::TwoTrackMinimumDistanceHelixHelix()
0011 : theH(nullptr), theG(nullptr), pointsUpdated(false), themaxjump(20), thesingjacI(1. / 0.1), themaxiter(4) {}
0012
0013 TwoTrackMinimumDistanceHelixHelix::~TwoTrackMinimumDistanceHelixHelix() {}
0014
0015 bool TwoTrackMinimumDistanceHelixHelix::updateCoeffs(const GlobalPoint& gpH, const GlobalPoint& gpG) {
0016 const double Bc2kH = theH->magneticField().inInverseGeV(gpH).z();
0017 const double Bc2kG = theG->magneticField().inInverseGeV(gpG).z();
0018 const double Ht = theH->momentum().perp();
0019 const double Gt = theG->momentum().perp();
0020
0021
0022 if (Ht == 0. || Gt == 0.) {
0023 edm::LogWarning("TwoTrackMinimumDistanceHelixHelix") << "transverse momentum of input trajectory is zero.";
0024 return true;
0025 };
0026
0027 if (theH->charge() == 0. || theG->charge() == 0.) {
0028 edm::LogWarning("TwoTrackMinimumDistanceHelixHelix") << "charge of input track is zero.";
0029 return true;
0030 };
0031
0032 if (Bc2kG == 0.) {
0033 edm::LogWarning("TwoTrackMinimumDistanceHelixHelix") << "magnetic field at point " << gpG << " is zero.";
0034 return true;
0035 };
0036
0037 if (Bc2kH == 0.) {
0038 edm::LogWarning("TwoTrackMinimumDistanceHelixHelix") << "magnetic field at point " << gpH << " is zero.";
0039 return true;
0040 };
0041
0042 theh = Ht / (theH->charge() * Bc2kH);
0043 thesinpH0 = -theH->momentum().y() / (Ht);
0044 thecospH0 = -theH->momentum().x() / (Ht);
0045 thetanlambdaH = -theH->momentum().z() / (Ht);
0046 thepH0 = atan2(thesinpH0, thecospH0);
0047
0048
0049
0050 theg = Gt / (theG->charge() * Bc2kG);
0051 thesinpG0 = -theG->momentum().y() / (Gt);
0052 thecospG0 = -theG->momentum().x() / (Gt);
0053 thetanlambdaG = -theG->momentum().z() / (Gt);
0054 thepG0 = atan2(thesinpG0, thecospG0);
0055
0056 thea = theH->position().x() - theG->position().x() + theg * thesinpG0 - theh * thesinpH0;
0057 theb = theH->position().y() - theG->position().y() - theg * thecospG0 + theh * thecospH0;
0058 thec1 = theh * thetanlambdaH * thetanlambdaH;
0059 thec2 = -theg * thetanlambdaG * thetanlambdaG;
0060 thed1 = -theg * thetanlambdaG * thetanlambdaH;
0061 thed2 = theh * thetanlambdaG * thetanlambdaH;
0062 thee1 = thetanlambdaH *
0063 (theH->position().z() - theG->position().z() - theh * thepH0 * thetanlambdaH + theg * thepG0 * thetanlambdaG);
0064 thee2 = thetanlambdaG *
0065 (theH->position().z() - theG->position().z() - theh * thepH0 * thetanlambdaH + theg * thepG0 * thetanlambdaG);
0066 return false;
0067 }
0068
0069 bool TwoTrackMinimumDistanceHelixHelix::oneIteration(double& dH, double& dG) {
0070 thesinpH = sin(thepH);
0071 thecospH = cos(thepH);
0072 thesinpG = sin(thepG);
0073 thecospG = cos(thepG);
0074
0075 const double A11 = theh * (-thesinpH * (thea - theg * thesinpG) + thecospH * (theb + theg * thecospG) + thec1);
0076 if (A11 < 0) {
0077 return true;
0078 };
0079 const double A22 = -theg * (-thesinpG * (thea + theh * thesinpH) + thecospG * (theb - theh * thecospH) + thec2);
0080 if (A22 < 0) {
0081 return true;
0082 };
0083 const double A12 = theh * (-theg * thecospG * thecospH - theg * thesinpH * thesinpG + thed1);
0084 const double A21 = -theg * (theh * thecospG * thecospH + theh * thesinpH * thesinpG + thed2);
0085 const double detaI = 1. / (A11 * A22 - A12 * A21);
0086 const double z1 = theh * (thecospH * (thea - theg * thesinpG) + thesinpH * (theb + theg * thecospG) + thec1 * thepH +
0087 thed1 * thepG + thee1);
0088 const double z2 = -theg * (thecospG * (thea + theh * thesinpH) + thesinpG * (theb - theh * thecospH) + thec2 * thepG +
0089 thed2 * thepH + thee2);
0090
0091 dH = (z1 * A22 - z2 * A12) * detaI;
0092 dG = (z2 * A11 - z1 * A21) * detaI;
0093 if (fabs(detaI) > thesingjacI) {
0094 return true;
0095 };
0096
0097 thepH -= dH;
0098 thepG -= dG;
0099 return false;
0100 }
0101
0102
0103
0104
0105
0106
0107
0108
0109
0110
0111
0112 bool TwoTrackMinimumDistanceHelixHelix::calculate(const GlobalTrajectoryParameters& G,
0113 const GlobalTrajectoryParameters& H,
0114 const float qual) {
0115 pointsUpdated = false;
0116 theG = &G;
0117 theH = &H;
0118 bool retval = false;
0119
0120 if (updateCoeffs(theG->position(), theH->position())) {
0121 finalPoints();
0122 return true;
0123 }
0124
0125 thepG = thepG0;
0126 thepH = thepH0;
0127
0128 int counter = 0;
0129 double pH = 0;
0130 double pG = 0;
0131 do {
0132 retval = oneIteration(pG, pH);
0133 if (edm::isNotFinite(pG) || edm::isNotFinite(pH))
0134 retval = true;
0135 if (counter++ > themaxiter)
0136 retval = true;
0137 } while ((!retval) && (fabs(pG) > qual || fabs(pH) > qual));
0138 if (fabs(theg * (thepG - thepG0)) > themaxjump)
0139 retval = true;
0140 if (fabs(theh * (thepH - thepH0)) > themaxjump)
0141 retval = true;
0142
0143 finalPoints();
0144
0145 return retval;
0146 }
0147
0148 void TwoTrackMinimumDistanceHelixHelix::finalPoints() {
0149 if (pointsUpdated)
0150 return;
0151 GlobalVector tmpG(sin(thepG) - thesinpG0, -cos(thepG) + thecospG0, thetanlambdaG * (thepG - thepG0));
0152 pointG = theG->position() + theg * tmpG;
0153 pathG = (thepG - thepG0) * (theg * sqrt(1 + thetanlambdaG * thetanlambdaG));
0154
0155 GlobalVector tmpH(sin(thepH) - thesinpH0, -cos(thepH) + thecospH0, thetanlambdaH * (thepH - thepH0));
0156 pointH = theH->position() + theh * tmpH;
0157 pathH = (thepH - thepH0) * (theh * sqrt(1 + thetanlambdaH * thetanlambdaH));
0158
0159 pointsUpdated = true;
0160 }
