PatternTools/src/TwoTrackMinimumDistanceHelixHelix.cc

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   //  thelambdaH=asin ( theH->momentum().z() / Hn );
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   // thelambdaG=asin ( theG->momentum().z()/( Gn ) );
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 bool TwoTrackMinimumDistanceHelixHelix::parallelTracks() const {
0104   return (fabs(theH->momentum().x() - theG->momentum().x()) < .00000001 )
0105     && (fabs(theH->momentum().y() - theG->momentum().y()) < .00000001 )
0106     && (fabs(theH->momentum().z() - theG->momentum().z()) < .00000001 )
0107     && (theH->charge()==theG->charge())
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 }