Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-06-24 02:11:23

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;
0037   Vec3D t2;
0038 
0039   Vec3D hn;
0040   Vec3D dx;
0041 
0042   for (int i = 0; i != 4; ++i) {
0043     t1[i] = p1.basicVector().v[i];
0044     t2[i] = p2.basicVector().v[i];
0045     hn[i] = hnf.basicVector().v[i];
0046     dx[i] = dxf.basicVector().v[i];
0047   }
0048 
0049   double gamma = dot(hn, t2);
0050   Vec3D an = cross3(hn, t2);
0051 
0052   typedef unsigned long long v4sl __attribute__((vector_size(32)));
0053   v4sl mask{1, 0, 5, 4};
0054   Vec4D t = __builtin_shuffle(t1, t2, mask);
0055   Vec4D tt = t * t;
0056 
0057   double au1 = std::sqrt(tt[0] + tt[1]);
0058   double au2 = std::sqrt(tt[2] + tt[3]);
0059   Vec4D au{-au1, au1, -au2, au2};
0060   Vec4D uu = t / au;
0061 
0062   Vec2D u1{uu[0], uu[1]};
0063   Vec2D u2{uu[2], uu[3]};
0064 
0065   Vec3D u13{uu[0], uu[1], 0, 0};
0066   Vec3D v1 = cross3(t1, u13);
0067   Vec3D u23{uu[2], uu[3], 0, 0};
0068   Vec3D v2 = cross3(t2, u23);
0069 
0070   // now prepare the transport matrix
0071   // pp only needed in high-p case (WA)
0072   //   double pp = 1./qbp;
0073   ////    double pp = fts.momentum().mag();
0074   // moved up (where -h.mag() is needed()
0075   //   double qp = q*pp;
0076 
0077   double anv = -dot(hn, u23);
0078   double anu = dot(hn, v2);
0079 
0080   double omcost = 1. - cost;
0081   double tmsint = theta - sint;
0082 
0083   Vec3D hu = cross3(hn, u13);
0084   Vec3D hv = cross3(hn, v1);
0085 
0086   //   1/p - doesn't change since |p1| = |p2|
0087   theJacobian(0, 0) = 1.;
0088   for (auto i = 1; i < 5; ++i)
0089     theJacobian(0, i) = 0.;
0090   //   lambda
0091 
0092   theJacobian(1, 0) = -qp * anv * dot(t2, dx);
0093 
0094   theJacobian(1, 1) = cost * dot(v1, v2) + sint * dot(hv, v2) + omcost * dot(hn, v1) * dot(hn, v2) +
0095                       anv * (-sint * dot(v1, t2) + omcost * dot(v1, an) - tmsint * gamma * dot(hn, v1));
0096 
0097   theJacobian(1, 2) = cost * dot2(u1, v2) + sint * dot(hu, v2) + omcost * dot2(hn, u1) * dot(hn, v2) +
0098                       anv * (-sint * dot2(u1, t2) + omcost * dot2(u1, an) - tmsint * gamma * dot2(hn, u1));
0099   theJacobian(1, 2) *= cosl0;
0100 
0101   theJacobian(1, 3) = -q * anv * dot2(u1, t2);
0102 
0103   theJacobian(1, 4) = -q * anv * dot(v1, t2);
0104 
0105   //   phi
0106 
0107   theJacobian(2, 0) = -qp * anu * cosl1 * dot(t2, dx);
0108 
0109   theJacobian(2, 1) = cost * dot(xy(v1), u2) + sint * dot(xy(hv), u2) + omcost * dot(hn, v1) * dot2(hn, u2) +
0110                       anu * (-sint * dot(v1, t2) + omcost * dot(v1, an) - tmsint * gamma * dot(hn, v1));
0111   theJacobian(2, 1) *= cosl1;
0112 
0113   theJacobian(2, 2) = cost * dot(u1, u2) + sint * dot2(hu, u2) + omcost * dot2(hn, u1) * dot2(hn, u2) +
0114                       anu * (-sint * dot2(u1, t2) + omcost * dot2(u1, an) - tmsint * gamma * dot2(hn, u1));
0115   theJacobian(2, 2) *= cosl1 * cosl0;
0116 
0117   theJacobian(2, 3) = -q * anu * cosl1 * dot2(u1, t2);
0118 
0119   theJacobian(2, 4) = -q * anu * cosl1 * dot(v1, t2);
0120 
0121   //   yt
0122 
0123   double overQ = 1. / q;
0124 
0125   theJacobian(3, 1) = (sint * dot(xy(v1), u2) + omcost * dot2(hv, u2) + tmsint * dot(hn, v1) * dot2(hn, u2)) * overQ;
0126 
0127   theJacobian(3, 2) =
0128       (sint * dot(u1, u2) + omcost * dot2(hu, u2) + tmsint * dot2(hn, u1) * dot2(hn, u2)) * (cosl0 * overQ);
0129 
0130   theJacobian(3, 3) = dot(u1, u2);
0131 
0132   theJacobian(3, 4) = dot(xy(v1), u2);
0133 
0134   //   zt
0135 
0136   theJacobian(4, 1) = (sint * dot(v1, v2) + omcost * dot(hv, v2) + tmsint * dot(hn, v1) * dot(hn, v2)) * overQ;
0137 
0138   theJacobian(4, 2) =
0139       (sint * dot(u1, xy(v2)) + omcost * dot(hu, v2) + tmsint * dot2(hn, u1) * dot(hn, v2)) * (cosl0 * overQ);
0140 
0141   theJacobian(4, 3) = dot2(u1, v2);
0142 
0143   theJacobian(4, 4) = dot(v1, v2);
0144 
0145   //double cutCriterion = abs(s/fts.momentum().mag());
0146   //  double cutCriterion = fabs(s/globalParameters.momentum().mag());
0147   double cutCriterion = std::abs(s * qbp);
0148   const double limit = 5.;  // valid for propagations with effectively float precision
0149 
0150   if (cutCriterion > limit) {
0151     double pp = 1. / qbp;
0152     theJacobian(3, 0) = pp * dot2(u2, dx);
0153     //    theJacobian(4,0) = -pp*dot(v2,dx);  // was wrong sign???
0154     /*
0155 It seems so.
0156 The effect was noticed analysing the distribution of reduced chi2 of general tracks
0157 vs eta. A +3% difference was noticed whem using the - sign instead of the plus,
0158 in the region .75 < |eta| < 1.5, in particular for <1GeV tracks.
0159 Indeed, the reduced chi2 is only one of the symptoms: many other quantities 
0160 were affected (momentum for example). In addition additional tracks were reconstructed.
0161 */
0162     theJacobian(4, 0) = pp * dot(v2, dx);
0163 
0164   } else {
0165     Vec3D hp1 = cross3(hn, t1);
0166     double temp1 = dot2(hp1, u2);
0167     Vec3D ghnmp = gamma * hn - t1;
0168     double temp2 = dot(xy(ghnmp), u2);
0169 
0170     double qps = qp * s;
0171     double h2 = qps * qbp;
0172     double h3 = (-1. / 8.) * h2;
0173 
0174     double secondOrder41 = 0.5 * temp1;
0175     double thirdOrder41 = (1. / 3.) * temp2;
0176     double fourthOrder41 = h3 * temp1;
0177     theJacobian(3, 0) = (s * qps) * (secondOrder41 + h2 * (thirdOrder41 + fourthOrder41));
0178 
0179     double temp3 = dot(hp1, v2);
0180     double temp4 = dot(ghnmp, v2);
0181 
0182     double secondOrder51 = 0.5 * temp3;
0183     double thirdOrder51 = (1. / 3.) * temp4;
0184     double fourthOrder51 = h3 * temp3;
0185     theJacobian(4, 0) = (s * qps) * (secondOrder51 + h2 * (thirdOrder51 + fourthOrder51));
0186   }
0187 
0188   // end of TRPRFN
0189 }