AnalyticalJacobians/src/AnalyticalCurvilinearJacobianEXT.icc

0001 // NOLINTNEXTLINE(misc-definitions-in-headers)
0002 void AnalyticalCurvilinearJacobian::computeFullJacobian(const GlobalTrajectoryParameters& globalParameters,
0003                                                         const GlobalPoint& x,
0004                                                         const GlobalVector& p,
0005                                                         const GlobalVector& h,
0006                                                         const double& s) {
0007   //GlobalVector p1 = fts.momentum().unit();
0008   GlobalVector p1 = globalParameters.momentum().unit();
0009   GlobalVector p2 = p.unit();
0010   //GlobalPoint xStart = fts.position();
0011   GlobalPoint xStart = globalParameters.position();
0012   GlobalVector dxf = xStart - x;
0013   //GlobalVector h  = MagneticField::inInverseGeV(xStart);
0014   // Martijn: field is now given as parameter.. GlobalVector h  = globalParameters.magneticFieldInInverseGeV(xStart);
0015
0016   //double qbp = fts.signedInverseMomentum();
0017   double qbp = globalParameters.signedInverseMomentum();
0018   double absS = s;
0019
0020   // calculate transport matrix
0021   // Origin: TRPRFN
0022   double cosl0 = p1.perp();
0023   double cosl1 = 1. / p2.perp();
0024
0025   // define average magnetic field and gradient
0026   // at initial point - inlike TRPRFN
0027   GlobalVector hnf = h.unit();
0028   double qp = -h.mag();
0029   //   double q = -h.mag()*qbp;
0030   double q = qp * qbp;
0031
0032   double theta = q * absS;
0033   double sint, cost;
0034   vdt::fast_sincos(theta, sint, cost);
0035
0036   Vec3D t1 = convert<Vec3D>(p1.basicVector().v);
0037   Vec3D t2 = convert<Vec3D>(p2.basicVector().v);
0038
0039   Vec3D hn = convert<Vec3D>(hnf.basicVector().v);
0040   Vec3D dx = convert<Vec3D>(dxf.basicVector().v);
0041
0042   double gamma = dot(hn, t2);
0043   Vec3D an = cross3(hn, t2);
0044
0045 #ifdef __clang__
0046   Vec4D t = __builtin_shufflevector(t1, t2, 1, 0, 5, 4);
0047 #else
0048   typedef unsigned long long v4sl __attribute__((vector_size(32)));
0049   v4sl mask{1, 0, 5, 4};
0050   Vec4D t = __builtin_shuffle(t1, t2, mask);
0051 #endif
0052   Vec4D tt = t * t;
0053
0054   Vec4D au{tt[1], tt[0], tt[3], tt[2]};
0055   au += tt;
0056   Vec4D neg{-1., 1., -1., 1.};
0057   for (int i = 0; i < 4; ++i)
0058     au[i] = neg[i] * std::sqrt(au[i]);
0059   /* above equivalent to this
0060   double au1 = std::sqrt(tt[0] + tt[1]);
0061   double au2 = std::sqrt(tt[2] + tt[3]);
0062   Vec4D au{-au1, au1, -au2, au2};
0063   */
0064   Vec4D uu = t / au;
0065
0066   Vec2D u1{uu[0], uu[1]};
0067   Vec2D u2{uu[2], uu[3]};
0068
0069   Vec3D u13{uu[0], uu[1], 0, 0};
0070   Vec3D v1 = cross3(t1, u13);
0071   Vec3D u23{uu[2], uu[3], 0, 0};
0072   Vec3D v2 = cross3(t2, u23);
0073
0074   // now prepare the transport matrix
0075   // pp only needed in high-p case (WA)
0076   //   double pp = 1./qbp;
0077   ////    double pp = fts.momentum().mag();
0078   // moved up (where -h.mag() is needed()
0079   //   double qp = q*pp;
0080
0081   double anv = -dot(hn, u23);
0082   double anu = dot(hn, v2);
0083
0084   double omcost = 1. - cost;
0085   double tmsint = theta - sint;
0086
0087   Vec3D hu = cross3(hn, u13);
0088   Vec3D hv = cross3(hn, v1);
0089
0090   //   1/p - doesn't change since |p1| = |p2|
0091   theJacobian(0, 0) = 1.;
0092   for (auto i = 1; i < 5; ++i)
0093     theJacobian(0, i) = 0.;
0094   //   lambda
0095
0096   theJacobian(1, 0) = -qp * anv * dot(t2, dx);
0097
0098   theJacobian(1, 1) = cost * dot(v1, v2) + sint * dot(hv, v2) + omcost * dot(hn, v1) * dot(hn, v2) +
0099                       anv * (-sint * dot(v1, t2) + omcost * dot(v1, an) - tmsint * gamma * dot(hn, v1));
0100
0101   theJacobian(1, 2) = cost * dot2(u1, v2) + sint * dot(hu, v2) + omcost * dot2(hn, u1) * dot(hn, v2) +
0102                       anv * (-sint * dot2(u1, t2) + omcost * dot2(u1, an) - tmsint * gamma * dot2(hn, u1));
0103   theJacobian(1, 2) *= cosl0;
0104
0105   theJacobian(1, 3) = -q * anv * dot2(u1, t2);
0106
0107   theJacobian(1, 4) = -q * anv * dot(v1, t2);
0108
0109   //   phi
0110
0111   theJacobian(2, 0) = -qp * anu * cosl1 * dot(t2, dx);
0112
0113   theJacobian(2, 1) = cost * dot(xy(v1), u2) + sint * dot(xy(hv), u2) + omcost * dot(hn, v1) * dot2(hn, u2) +
0114                       anu * (-sint * dot(v1, t2) + omcost * dot(v1, an) - tmsint * gamma * dot(hn, v1));
0115   theJacobian(2, 1) *= cosl1;
0116
0117   theJacobian(2, 2) = cost * dot(u1, u2) + sint * dot2(hu, u2) + omcost * dot2(hn, u1) * dot2(hn, u2) +
0118                       anu * (-sint * dot2(u1, t2) + omcost * dot2(u1, an) - tmsint * gamma * dot2(hn, u1));
0119   theJacobian(2, 2) *= cosl1 * cosl0;
0120
0121   theJacobian(2, 3) = -q * anu * cosl1 * dot2(u1, t2);
0122
0123   theJacobian(2, 4) = -q * anu * cosl1 * dot(v1, t2);
0124
0125   //   yt
0126
0127   double overQ = 1. / q;
0128
0129   theJacobian(3, 1) = (sint * dot(xy(v1), u2) + omcost * dot2(hv, u2) + tmsint * dot(hn, v1) * dot2(hn, u2)) * overQ;
0130
0131   theJacobian(3, 2) =
0132       (sint * dot(u1, u2) + omcost * dot2(hu, u2) + tmsint * dot2(hn, u1) * dot2(hn, u2)) * (cosl0 * overQ);
0133
0134   theJacobian(3, 3) = dot(u1, u2);
0135
0136   theJacobian(3, 4) = dot(xy(v1), u2);
0137
0138   //   zt
0139
0140   theJacobian(4, 1) = (sint * dot(v1, v2) + omcost * dot(hv, v2) + tmsint * dot(hn, v1) * dot(hn, v2)) * overQ;
0141
0142   theJacobian(4, 2) =
0143       (sint * dot(u1, xy(v2)) + omcost * dot(hu, v2) + tmsint * dot2(hn, u1) * dot(hn, v2)) * (cosl0 * overQ);
0144
0145   theJacobian(4, 3) = dot2(u1, v2);
0146
0147   theJacobian(4, 4) = dot(v1, v2);
0148
0149   //double cutCriterion = abs(s/fts.momentum().mag());
0150   //  double cutCriterion = fabs(s/globalParameters.momentum().mag());
0151   double cutCriterion = std::abs(s * qbp);
0152   const double limit = 5.;  // valid for propagations with effectively float precision
0153
0154   if (cutCriterion > limit) {
0155     double pp = 1. / qbp;
0156     theJacobian(3, 0) = pp * dot2(u2, dx);
0157     //    theJacobian(4,0) = -pp*dot(v2,dx);  // was wrong sign???
0158     /*
0159 It seems so.
0160 The effect was noticed analysing the distribution of reduced chi2 of general tracks
0161 vs eta. A +3% difference was noticed whem using the - sign instead of the plus,
0162 in the region .75 < |eta| < 1.5, in particular for <1GeV tracks.
0163 Indeed, the reduced chi2 is only one of the symptoms: many other quantities
0164 were affected (momentum for example). In addition additional tracks were reconstructed.
0165 */
0166     theJacobian(4, 0) = pp * dot(v2, dx);
0167
0168   } else {
0169     Vec3D hp1 = cross3(hn, t1);
0170     double temp1 = dot2(hp1, u2);
0171     Vec3D ghnmp = gamma * hn - t1;
0172     double temp2 = dot(xy(ghnmp), u2);
0173
0174     double qps = qp * s;
0175     double h2 = qps * qbp;
0176     double h3 = (-1. / 8.) * h2;
0177
0178     double secondOrder41 = 0.5 * temp1;
0179     double thirdOrder41 = (1. / 3.) * temp2;
0180     double fourthOrder41 = h3 * temp1;
0181     theJacobian(3, 0) = (s * qps) * (secondOrder41 + h2 * (thirdOrder41 + fourthOrder41));
0182
0183     double temp3 = dot(hp1, v2);
0184     double temp4 = dot(ghnmp, v2);
0185
0186     double secondOrder51 = 0.5 * temp3;
0187     double thirdOrder51 = (1. / 3.) * temp4;
0188     double fourthOrder51 = h3 * temp3;
0189     theJacobian(4, 0) = (s * qps) * (secondOrder51 + h2 * (thirdOrder51 + fourthOrder51));
0190   }
0191
0192   // end of TRPRFN
0193 }