MkFitCore/src/PropagationMPlexEndcap.cc

0001 #include "FWCore/Utilities/interface/CMSUnrollLoop.h"
0002 #include "RecoTracker/MkFitCore/interface/PropagationConfig.h"
0003 #include "RecoTracker/MkFitCore/interface/Config.h"
0004 #include "RecoTracker/MkFitCore/interface/TrackerInfo.h"
0005
0006 #include "PropagationMPlex.h"
0007
0008 //#define DEBUG
0009 #include "Debug.h"
0010
0011 namespace {
0012   using namespace mkfit;
0013
0014   void MultHelixPropEndcap(const MPlexLL& A, const MPlexLS& B, MPlexLL& C) {
0015     // C = A * B
0016
0017     typedef float T;
0018     const Matriplex::idx_t N = NN;
0019
0020     const T* a = A.fArray;
0021     ASSUME_ALIGNED(a, 64);
0022     const T* b = B.fArray;
0023     ASSUME_ALIGNED(b, 64);
0024     T* c = C.fArray;
0025     ASSUME_ALIGNED(c, 64);
0026
0027 #include "MultHelixPropEndcap.ah"
0028   }
0029
0030   void MultHelixPropTranspEndcap(const MPlexLL& A, const MPlexLL& B, MPlexLS& C) {
0031     // C = B * AT;
0032
0033     typedef float T;
0034     const Matriplex::idx_t N = NN;
0035
0036     const T* a = A.fArray;
0037     ASSUME_ALIGNED(a, 64);
0038     const T* b = B.fArray;
0039     ASSUME_ALIGNED(b, 64);
0040     T* c = C.fArray;
0041     ASSUME_ALIGNED(c, 64);
0042
0043 #include "MultHelixPropTranspEndcap.ah"
0044   }
0045
0046 }  // namespace
0047
0048 // ============================================================================
0049 // BEGIN STUFF FROM PropagationMPlex.icc
0050 namespace {}
0051
0052 // END STUFF FROM PropagationMPlex.icc
0053 // ============================================================================
0054
0055 namespace mkfit {
0056
0057   void propagateHelixToZMPlex(const MPlexLS& inErr,
0058                               const MPlexLV& inPar,
0059                               const MPlexQI& inChg,
0060                               const MPlexQF& msZ,
0061                               MPlexLS& outErr,
0062                               MPlexLV& outPar,
0063                               MPlexQI& outFailFlag,
0064                               const int N_proc,
0065                               const PropagationFlags& pflags,
0066                               const MPlexQI* noMatEffPtr) {
0067     // debug = true;
0068
0069     outErr = inErr;
0070     outPar = inPar;
0071
0072     MPlexLL errorProp;
0073
0074     //helixAtZ_new(inPar, inChg, msZ, outPar, errorProp, outFailFlag, N_proc, pflags);
0075     helixAtZ(inPar, inChg, msZ, outPar, errorProp, outFailFlag, N_proc, pflags);
0076
0077 #ifdef DEBUG
0078     if (debug && g_debug) {
0079       for (int kk = 0; kk < N_proc; ++kk) {
0080         dprintf("inPar %d\n", kk);
0081         for (int i = 0; i < 6; ++i) {
0082           dprintf("%8f ", inPar.constAt(kk, i, 0));
0083         }
0084         dprintf("\n");
0085
0086         dprintf("inErr %d\n", kk);
0087         for (int i = 0; i < 6; ++i) {
0088           for (int j = 0; j < 6; ++j)
0089             dprintf("%8f ", inErr.constAt(kk, i, j));
0090           dprintf("\n");
0091         }
0092         dprintf("\n");
0093
0094         dprintf("errorProp %d\n", kk);
0095         for (int i = 0; i < 6; ++i) {
0096           for (int j = 0; j < 6; ++j)
0097             dprintf("%8f ", errorProp.At(kk, i, j));
0098           dprintf("\n");
0099         }
0100         dprintf("\n");
0101       }
0102     }
0103 #endif
0104
0105 #ifdef DEBUG
0106     if (debug && g_debug) {
0107       for (int kk = 0; kk < N_proc; ++kk) {
0108         dprintf("outErr %d\n", kk);
0109         for (int i = 0; i < 6; ++i) {
0110           for (int j = 0; j < 6; ++j)
0111             dprintf("%8f ", outErr.constAt(kk, i, j));
0112           dprintf("\n");
0113         }
0114         dprintf("\n");
0115       }
0116     }
0117 #endif
0118
0119     // Matriplex version of: result.errors = ROOT::Math::Similarity(errorProp, outErr);
0120     MPlexLL temp;
0121     MultHelixPropEndcap(errorProp, outErr, temp);
0122     MultHelixPropTranspEndcap(errorProp, temp, outErr);
0123     // can replace with: MultHelixPropFull(errorProp, outErr, temp); MultHelixPropTranspFull(errorProp, temp, outErr);
0124
0125     if (pflags.apply_material) {
0126       MPlexQF hitsRl;
0127       MPlexQF hitsXi;
0128       MPlexQF propSign;
0129
0130       const TrackerInfo& tinfo = *pflags.tracker_info;
0131
0132 #pragma omp simd
0133       for (int n = 0; n < NN; ++n) {
0134         if (n >= N_proc || (noMatEffPtr && noMatEffPtr->constAt(n, 0, 0))) {
0135           hitsRl(n, 0, 0) = 0.f;
0136           hitsXi(n, 0, 0) = 0.f;
0137         } else {
0138           const float hypo = std::hypot(outPar(n, 0, 0), outPar(n, 1, 0));
0139           auto mat = tinfo.material_checked(std::abs(msZ(n, 0, 0)), hypo);
0140           hitsRl(n, 0, 0) = mat.radl;
0141           hitsXi(n, 0, 0) = mat.bbxi;
0142         }
0143         if (n < N_proc) {
0144           const float zout = msZ.constAt(n, 0, 0);
0145           const float zin = inPar.constAt(n, 2, 0);
0146           propSign(n, 0, 0) = (std::abs(zout) > std::abs(zin) ? 1.f : -1.f);
0147         }
0148       }
0149       MPlexHV plNrm;
0150 #pragma omp simd
0151       for (int n = 0; n < NN; ++n) {
0152         plNrm(n, 0, 0) = 0.f;
0153         plNrm(n, 1, 0) = 0.f;
0154         plNrm(n, 2, 0) = 1.f;
0155       }
0156       applyMaterialEffects(hitsRl, hitsXi, propSign, plNrm, outErr, outPar, N_proc);
0157 #ifdef DEBUG
0158       if (debug && g_debug) {
0159         for (int kk = 0; kk < N_proc; ++kk) {
0160           dprintf("propSign %d\n", kk);
0161           for (int i = 0; i < 1; ++i) {
0162             dprintf("%8f ", propSign.constAt(kk, i, 0));
0163           }
0164           dprintf("\n");
0165           dprintf("plNrm %d\n", kk);
0166           for (int i = 0; i < 3; ++i) {
0167             dprintf("%8f ", plNrm.constAt(kk, i, 0));
0168           }
0169           dprintf("\n");
0170           dprintf("outErr(after material) %d\n", kk);
0171           for (int i = 0; i < 6; ++i) {
0172             for (int j = 0; j < 6; ++j)
0173               dprintf("%8f ", outErr.constAt(kk, i, j));
0174             dprintf("\n");
0175           }
0176           dprintf("\n");
0177         }
0178       }
0179 #endif
0180     }
0181
0182     squashPhiMPlex(outPar, N_proc);  // ensure phi is between |pi|
0183
0184     // PROP-FAIL-ENABLE To keep physics changes minimal, we always restore the
0185     // state to input when propagation fails -- as was the default before.
0186     // if (pflags.copy_input_state_on_fail) {
0187     for (int i = 0; i < N_proc; ++i) {
0188       if (outFailFlag(i, 0, 0)) {
0189         outPar.copySlot(i, inPar);
0190         outErr.copySlot(i, inErr);
0191       }
0192     }
0193     // }
0194   }
0195
0196   void helixAtZ(const MPlexLV& inPar,
0197                 const MPlexQI& inChg,
0198                 const MPlexQF& msZ,
0199                 MPlexLV& outPar,
0200                 MPlexLL& errorProp,
0201                 MPlexQI& outFailFlag,
0202                 const int N_proc,
0203                 const PropagationFlags& pflags) {
0204     errorProp.setVal(0.f);
0205     outFailFlag.setVal(0.f);
0206
0207     // debug = true;
0208 #pragma omp simd
0209     for (int n = 0; n < NN; ++n) {
0210       //initialize erroProp to identity matrix, except element 2,2 which is zero
0211       errorProp(n, 0, 0) = 1.f;
0212       errorProp(n, 1, 1) = 1.f;
0213       errorProp(n, 3, 3) = 1.f;
0214       errorProp(n, 4, 4) = 1.f;
0215       errorProp(n, 5, 5) = 1.f;
0216     }
0217     float zout[NN];
0218     float zin[NN];
0219     float ipt[NN];
0220     float phiin[NN];
0221     float theta[NN];
0222 #pragma omp simd
0223     for (int n = 0; n < NN; ++n) {
0224       //initialize erroProp to identity matrix, except element 2,2 which is zero
0225       zout[n] = msZ.constAt(n, 0, 0);
0226       zin[n] = inPar.constAt(n, 2, 0);
0227       ipt[n] = inPar.constAt(n, 3, 0);
0228       phiin[n] = inPar.constAt(n, 4, 0);
0229       theta[n] = inPar.constAt(n, 5, 0);
0230     }
0231
0232     float k[NN];
0233     if (pflags.use_param_b_field) {
0234 #pragma omp simd
0235       for (int n = 0; n < NN; ++n) {
0236         k[n] = inChg.constAt(n, 0, 0) * 100.f /
0237                (-Const::sol * Config::bFieldFromZR(zin[n], hipo(inPar.constAt(n, 0, 0), inPar.constAt(n, 1, 0))));
0238       }
0239     } else {
0240 #pragma omp simd
0241       for (int n = 0; n < NN; ++n) {
0242         k[n] = inChg.constAt(n, 0, 0) * 100.f / (-Const::sol * Config::Bfield);
0243       }
0244     }
0245
0246     float kinv[NN];
0247 #pragma omp simd
0248     for (int n = 0; n < NN; ++n) {
0249       kinv[n] = 1.f / k[n];
0250     }
0251
0252 #pragma omp simd
0253     for (int n = 0; n < NN; ++n) {
0254       dprint_np(n,
0255                 std::endl
0256                     << "input parameters"
0257                     << " inPar.constAt(n, 0, 0)=" << std::setprecision(9) << inPar.constAt(n, 0, 0)
0258                     << " inPar.constAt(n, 1, 0)=" << std::setprecision(9) << inPar.constAt(n, 1, 0)
0259                     << " inPar.constAt(n, 2, 0)=" << std::setprecision(9) << inPar.constAt(n, 2, 0)
0260                     << " inPar.constAt(n, 3, 0)=" << std::setprecision(9) << inPar.constAt(n, 3, 0)
0261                     << " inPar.constAt(n, 4, 0)=" << std::setprecision(9) << inPar.constAt(n, 4, 0)
0262                     << " inPar.constAt(n, 5, 0)=" << std::setprecision(9) << inPar.constAt(n, 5, 0)
0263                     << " inChg.constAt(n, 0, 0)=" << std::setprecision(9) << inChg.constAt(n, 0, 0));
0264     }
0265 #pragma omp simd
0266     for (int n = 0; n < NN; ++n) {
0267       dprint_np(n,
0268                 "propagation start, dump parameters"
0269                     << std::endl
0270                     << "pos = " << inPar.constAt(n, 0, 0) << " " << inPar.constAt(n, 1, 0) << " "
0271                     << inPar.constAt(n, 2, 0) << std::endl
0272                     << "mom (cart) = " << std::cos(inPar.constAt(n, 4, 0)) / inPar.constAt(n, 3, 0) << " "
0273                     << std::sin(inPar.constAt(n, 4, 0)) / inPar.constAt(n, 3, 0) << " "
0274                     << 1. / (inPar.constAt(n, 3, 0) * tan(inPar.constAt(n, 5, 0))) << " r="
0275                     << std::sqrt(inPar.constAt(n, 0, 0) * inPar.constAt(n, 0, 0) +
0276                                  inPar.constAt(n, 1, 0) * inPar.constAt(n, 1, 0))
0277                     << " pT=" << 1. / std::abs(inPar.constAt(n, 3, 0)) << " q=" << inChg.constAt(n, 0, 0)
0278                     << " targetZ=" << msZ.constAt(n, 0, 0) << std::endl);
0279     }
0280
0281     float pt[NN];
0282 #pragma omp simd
0283     for (int n = 0; n < NN; ++n) {
0284       pt[n] = 1.f / ipt[n];
0285     }
0286
0287     //no trig approx here, phi can be large
0288     float cosP[NN];
0289     float sinP[NN];
0290 #pragma omp simd
0291     for (int n = 0; n < NN; ++n) {
0292       cosP[n] = std::cos(phiin[n]);
0293     }
0294
0295 #pragma omp simd
0296     for (int n = 0; n < NN; ++n) {
0297       sinP[n] = std::sin(phiin[n]);
0298     }
0299
0300     float cosT[NN];
0301     float sinT[NN];
0302 #pragma omp simd
0303     for (int n = 0; n < NN; ++n) {
0304       cosT[n] = std::cos(theta[n]);
0305     }
0306
0307 #pragma omp simd
0308     for (int n = 0; n < NN; ++n) {
0309       sinT[n] = std::sin(theta[n]);
0310     }
0311
0312     float tanT[NN];
0313     float icos2T[NN];
0314     float pxin[NN];
0315     float pyin[NN];
0316 #pragma omp simd
0317     for (int n = 0; n < NN; ++n) {
0318       tanT[n] = sinT[n] / cosT[n];
0319       icos2T[n] = 1.f / (cosT[n] * cosT[n]);
0320       pxin[n] = cosP[n] * pt[n];
0321       pyin[n] = sinP[n] * pt[n];
0322     }
0323
0324     float deltaZ[NN];
0325     float alpha[NN];
0326 #pragma omp simd
0327     for (int n = 0; n < NN; ++n) {
0328       deltaZ[n] = zout[n] - zin[n];
0329       alpha[n] = deltaZ[n] * tanT[n] * ipt[n] * kinv[n];
0330     }
0331
0332     float cosahTmp[NN];
0333     float sinahTmp[NN];
0334     if constexpr (Config::useTrigApprox) {
0335 #if !defined(__INTEL_COMPILER)
0336 #pragma omp simd
0337 #endif
0338       for (int n = 0; n < NN; ++n) {
0339         sincos4(alpha[n] * 0.5f, sinahTmp[n], cosahTmp[n]);
0340       }
0341     } else {
0342 #if !defined(__INTEL_COMPILER)
0343 #pragma omp simd
0344 #endif
0345       for (int n = 0; n < NN; ++n) {
0346         cosahTmp[n] = std::cos(alpha[n] * 0.5f);
0347       }
0348 #if !defined(__INTEL_COMPILER)
0349 #pragma omp simd
0350 #endif
0351       for (int n = 0; n < NN; ++n) {
0352         sinahTmp[n] = std::sin(alpha[n] * 0.5f);
0353       }
0354     }
0355
0356     float cosah[NN];
0357     float sinah[NN];
0358     float cosa[NN];
0359     float sina[NN];
0360 #pragma omp simd
0361     for (int n = 0; n < NN; ++n) {
0362       cosah[n] = cosahTmp[n];
0363       sinah[n] = sinahTmp[n];
0364       cosa[n] = 1.f - 2.f * sinah[n] * sinah[n];
0365       sina[n] = 2.f * sinah[n] * cosah[n];
0366     }
0367
0368 //update parameters
0369 #pragma omp simd
0370     for (int n = 0; n < NN; ++n) {
0371       outPar.At(n, 0, 0) = outPar.At(n, 0, 0) + 2.f * k[n] * sinah[n] * (pxin[n] * cosah[n] - pyin[n] * sinah[n]);
0372       outPar.At(n, 1, 0) = outPar.At(n, 1, 0) + 2.f * k[n] * sinah[n] * (pyin[n] * cosah[n] + pxin[n] * sinah[n]);
0373       outPar.At(n, 2, 0) = zout[n];
0374       outPar.At(n, 4, 0) = phiin[n] + alpha[n];
0375     }
0376
0377 #pragma omp simd
0378     for (int n = 0; n < NN; ++n) {
0379       dprint_np(n,
0380                 "propagation to Z end (OLD), dump parameters\n"
0381                     << "   pos = " << outPar(n, 0, 0) << " " << outPar(n, 1, 0) << " " << outPar(n, 2, 0) << "\t\t r="
0382                     << std::sqrt(outPar(n, 0, 0) * outPar(n, 0, 0) + outPar(n, 1, 0) * outPar(n, 1, 0)) << std::endl
0383                     << "   mom = " << outPar(n, 3, 0) << " " << outPar(n, 4, 0) << " " << outPar(n, 5, 0) << std::endl
0384                     << " cart= " << std::cos(outPar(n, 4, 0)) / outPar(n, 3, 0) << " "
0385                     << std::sin(outPar(n, 4, 0)) / outPar(n, 3, 0) << " "
0386                     << 1. / (outPar(n, 3, 0) * tan(outPar(n, 5, 0))) << "\t\tpT=" << 1. / std::abs(outPar(n, 3, 0))
0387                     << std::endl);
0388     }
0389
0390     float pxcaMpysa[NN];
0391 #pragma omp simd
0392     for (int n = 0; n < NN; ++n) {
0393       pxcaMpysa[n] = pxin[n] * cosa[n] - pyin[n] * sina[n];
0394     }
0395
0396 #pragma omp simd
0397     for (int n = 0; n < NN; ++n) {
0398       errorProp(n, 0, 2) = -tanT[n] * ipt[n] * pxcaMpysa[n];
0399       errorProp(n, 0, 3) =
0400           k[n] * pt[n] * pt[n] *
0401           (cosP[n] * (alpha[n] * cosa[n] - sina[n]) + sinP[n] * 2.f * sinah[n] * (sinah[n] - alpha[n] * cosah[n]));
0402       errorProp(n, 0, 4) = -2.f * k[n] * pt[n] * sinah[n] * (sinP[n] * cosah[n] + cosP[n] * sinah[n]);
0403       errorProp(n, 0, 5) = deltaZ[n] * ipt[n] * pxcaMpysa[n] * icos2T[n];
0404     }
0405
0406     float pycaPpxsa[NN];
0407 #pragma omp simd
0408     for (int n = 0; n < NN; ++n) {
0409       pycaPpxsa[n] = pyin[n] * cosa[n] + pxin[n] * sina[n];
0410     }
0411
0412 #pragma omp simd
0413     for (int n = 0; n < NN; ++n) {
0414       errorProp(n, 1, 2) = -tanT[n] * ipt[n] * pycaPpxsa[n];
0415       errorProp(n, 1, 3) =
0416           k[n] * pt[n] * pt[n] *
0417           (sinP[n] * (alpha[n] * cosa[n] - sina[n]) - cosP[n] * 2.f * sinah[n] * (sinah[n] - alpha[n] * cosah[n]));
0418       errorProp(n, 1, 4) = 2.f * k[n] * pt[n] * sinah[n] * (cosP[n] * cosah[n] - sinP[n] * sinah[n]);
0419       errorProp(n, 1, 5) = deltaZ[n] * ipt[n] * pycaPpxsa[n] * icos2T[n];
0420     }
0421
0422 #pragma omp simd
0423     for (int n = 0; n < NN; ++n) {
0424       errorProp(n, 4, 2) = -ipt[n] * tanT[n] * kinv[n];
0425       errorProp(n, 4, 3) = tanT[n] * deltaZ[n] * kinv[n];
0426       errorProp(n, 4, 5) = ipt[n] * deltaZ[n] * kinv[n] * icos2T[n];
0427     }
0428
0429 #pragma omp simd
0430     for (int n = 0; n < NN; ++n) {
0431       dprint_np(
0432           n,
0433           "propagation end, dump parameters"
0434               << std::endl
0435               << "pos = " << outPar.At(n, 0, 0) << " " << outPar.At(n, 1, 0) << " " << outPar.At(n, 2, 0) << std::endl
0436               << "mom (cart) = " << std::cos(outPar.At(n, 4, 0)) / outPar.At(n, 3, 0) << " "
0437               << std::sin(outPar.At(n, 4, 0)) / outPar.At(n, 3, 0) << " "
0438               << 1. / (outPar.At(n, 3, 0) * tan(outPar.At(n, 5, 0)))
0439               << " r=" << std::sqrt(outPar.At(n, 0, 0) * outPar.At(n, 0, 0) + outPar.At(n, 1, 0) * outPar.At(n, 1, 0))
0440               << " pT=" << 1. / std::abs(outPar.At(n, 3, 0)) << std::endl);
0441     }
0442
0443     // PROP-FAIL-ENABLE Disabled to keep physics changes minimal.
0444     // To be reviewed, enabled and processed accordingly elsewhere.
0445     /*
0446     // Check for errors, set fail-flag.
0447     for (int n = 0; n < NN; ++n) {
0448       // We propagate for alpha: mark fail when prop angle more than pi/2
0449       if (std::abs(alpha[n]) > 1.57) {
0450         dprintf("helixAtZ: more than quarter turn, alpha = %f\n", alpha[n]);
0451         outFailFlag[n] = 1;
0452       } else {
0453         // Have we reached desired z? We can't know, we copy desired z to actual z.
0454         // Are we close to apex? Same condition as in propToR, 12.5 deg, cos(78.5deg) = 0.2
0455         float dotp = (outPar.At(n, 0, 0) * std::cos(outPar.At(n, 4, 0)) +
0456                       outPar.At(n, 1, 0) * std::sin(outPar.At(n, 4, 0))) /
0457                      std::hypot(outPar.At(n, 0, 0), outPar.At(n, 1, 0));
0458         if (dotp < 0.2 || dotp < 0) {
0459           dprintf("helixAtZ: dot product bad, dotp = %f\n", dotp);
0460           outFailFlag[n] = 1;
0461         }
0462       }
0463     }
0464     */
0465
0466 #ifdef DEBUG
0467     if (debug && g_debug) {
0468       for (int n = 0; n < N_proc; ++n) {
0469         dmutex_guard;
0470         std::cout << n << ": jacobian" << std::endl;
0471         printf("%5f %5f %5f %5f %5f %5f\n",
0472                errorProp(n, 0, 0),
0473                errorProp(n, 0, 1),
0474                errorProp(n, 0, 2),
0475                errorProp(n, 0, 3),
0476                errorProp(n, 0, 4),
0477                errorProp(n, 0, 5));
0478         printf("%5f %5f %5f %5f %5f %5f\n",
0479                errorProp(n, 1, 0),
0480                errorProp(n, 1, 1),
0481                errorProp(n, 1, 2),
0482                errorProp(n, 1, 3),
0483                errorProp(n, 1, 4),
0484                errorProp(n, 1, 5));
0485         printf("%5f %5f %5f %5f %5f %5f\n",
0486                errorProp(n, 2, 0),
0487                errorProp(n, 2, 1),
0488                errorProp(n, 2, 2),
0489                errorProp(n, 2, 3),
0490                errorProp(n, 2, 4),
0491                errorProp(n, 2, 5));
0492         printf("%5f %5f %5f %5f %5f %5f\n",
0493                errorProp(n, 3, 0),
0494                errorProp(n, 3, 1),
0495                errorProp(n, 3, 2),
0496                errorProp(n, 3, 3),
0497                errorProp(n, 3, 4),
0498                errorProp(n, 3, 5));
0499         printf("%5f %5f %5f %5f %5f %5f\n",
0500                errorProp(n, 4, 0),
0501                errorProp(n, 4, 1),
0502                errorProp(n, 4, 2),
0503                errorProp(n, 4, 3),
0504                errorProp(n, 4, 4),
0505                errorProp(n, 4, 5));
0506         printf("%5f %5f %5f %5f %5f %5f\n",
0507                errorProp(n, 5, 0),
0508                errorProp(n, 5, 1),
0509                errorProp(n, 5, 2),
0510                errorProp(n, 5, 3),
0511                errorProp(n, 5, 4),
0512                errorProp(n, 5, 5));
0513       }
0514     }
0515 #endif
0516   }
0517
0518 }  // namespace mkfit