Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-01-18 03:42:21 

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 #if !defined(__clang__)
0133 #pragma omp simd
0134 #endif
0135       for (int n = 0; n < NN; ++n) {
0136         if (n >= N_proc || (noMatEffPtr && noMatEffPtr->constAt(n, 0, 0))) {
0137           hitsRl(n, 0, 0) = 0.f;
0138           hitsXi(n, 0, 0) = 0.f;
0139         } else {
0140           const float hypo = hipo(outPar(n, 0, 0), outPar(n, 1, 0));
0141           const auto mat = tinfo.material_checked(std::abs(msZ(n, 0, 0)), hypo);
0142           hitsRl(n, 0, 0) = mat.radl;
0143           hitsXi(n, 0, 0) = mat.bbxi;
0144         }
0145         if (n < N_proc) {
0146           const float zout = msZ.constAt(n, 0, 0);
0147           const float zin = inPar.constAt(n, 2, 0);
0148           propSign(n, 0, 0) = (std::abs(zout) > std::abs(zin) ? 1.f : -1.f);
0149         }
0150       }
0151       MPlexHV plNrm;
0152 #pragma omp simd
0153       for (int n = 0; n < NN; ++n) {
0154         plNrm(n, 0, 0) = 0.f;
0155         plNrm(n, 1, 0) = 0.f;
0156         plNrm(n, 2, 0) = 1.f;
0157       }
0158       applyMaterialEffects(hitsRl, hitsXi, propSign, plNrm, outErr, outPar, N_proc);
0159 #ifdef DEBUG
0160       if (debug && g_debug) {
0161         for (int kk = 0; kk < N_proc; ++kk) {
0162           dprintf("propSign %d\n", kk);
0163           for (int i = 0; i < 1; ++i) {
0164             dprintf("%8f ", propSign.constAt(kk, i, 0));
0165           }
0166           dprintf("\n");
0167           dprintf("plNrm %d\n", kk);
0168           for (int i = 0; i < 3; ++i) {
0169             dprintf("%8f ", plNrm.constAt(kk, i, 0));
0170           }
0171           dprintf("\n");
0172           dprintf("outErr(after material) %d\n", kk);
0173           for (int i = 0; i < 6; ++i) {
0174             for (int j = 0; j < 6; ++j)
0175               dprintf("%8f ", outErr.constAt(kk, i, j));
0176             dprintf("\n");
0177           }
0178           dprintf("\n");
0179         }
0180       }
0181 #endif
0182     }
0183 
0184     squashPhiMPlex(outPar, N_proc);  // ensure phi is between |pi|
0185 
0186     // PROP-FAIL-ENABLE To keep physics changes minimal, we always restore the
0187     // state to input when propagation fails -- as was the default before.
0188     // if (pflags.copy_input_state_on_fail) {
0189     for (int i = 0; i < N_proc; ++i) {
0190       if (outFailFlag(i, 0, 0)) {
0191         outPar.copySlot(i, inPar);
0192         outErr.copySlot(i, inErr);
0193       }
0194     }
0195     // }
0196   }
0197 
0198   void helixAtZ(const MPlexLV& inPar,
0199                 const MPlexQI& inChg,
0200                 const MPlexQF& msZ,
0201                 MPlexLV& outPar,
0202                 MPlexLL& errorProp,
0203                 MPlexQI& outFailFlag,
0204                 const int N_proc,
0205                 const PropagationFlags& pflags) {
0206     errorProp.setVal(0.f);
0207     outFailFlag.setVal(0.f);
0208 
0209     // debug = true;
0210 #pragma omp simd
0211     for (int n = 0; n < NN; ++n) {
0212       //initialize erroProp to identity matrix, except element 2,2 which is zero
0213       errorProp(n, 0, 0) = 1.f;
0214       errorProp(n, 1, 1) = 1.f;
0215       errorProp(n, 3, 3) = 1.f;
0216       errorProp(n, 4, 4) = 1.f;
0217       errorProp(n, 5, 5) = 1.f;
0218     }
0219     float zout[NN];
0220     float zin[NN];
0221     float ipt[NN];
0222     float phiin[NN];
0223     float theta[NN];
0224 #pragma omp simd
0225     for (int n = 0; n < NN; ++n) {
0226       //initialize erroProp to identity matrix, except element 2,2 which is zero
0227       zout[n] = msZ.constAt(n, 0, 0);
0228       zin[n] = inPar.constAt(n, 2, 0);
0229       ipt[n] = inPar.constAt(n, 3, 0);
0230       phiin[n] = inPar.constAt(n, 4, 0);
0231       theta[n] = inPar.constAt(n, 5, 0);
0232     }
0233 
0234     float k[NN];
0235     if (pflags.use_param_b_field) {
0236 #pragma omp simd
0237       for (int n = 0; n < NN; ++n) {
0238         k[n] = inChg.constAt(n, 0, 0) * 100.f /
0239                (-Const::sol * Config::bFieldFromZR(zin[n], hipo(inPar.constAt(n, 0, 0), inPar.constAt(n, 1, 0))));
0240       }
0241     } else {
0242 #pragma omp simd
0243       for (int n = 0; n < NN; ++n) {
0244         k[n] = inChg.constAt(n, 0, 0) * 100.f / (-Const::sol * Config::Bfield);
0245       }
0246     }
0247 
0248     float kinv[NN];
0249 #pragma omp simd
0250     for (int n = 0; n < NN; ++n) {
0251       kinv[n] = 1.f / k[n];
0252     }
0253 
0254 #pragma omp simd
0255     for (int n = 0; n < NN; ++n) {
0256       dprint_np(n,
0257                 std::endl
0258                     << "input parameters"
0259                     << " inPar.constAt(n, 0, 0)=" << std::setprecision(9) << inPar.constAt(n, 0, 0)
0260                     << " inPar.constAt(n, 1, 0)=" << std::setprecision(9) << inPar.constAt(n, 1, 0)
0261                     << " inPar.constAt(n, 2, 0)=" << std::setprecision(9) << inPar.constAt(n, 2, 0)
0262                     << " inPar.constAt(n, 3, 0)=" << std::setprecision(9) << inPar.constAt(n, 3, 0)
0263                     << " inPar.constAt(n, 4, 0)=" << std::setprecision(9) << inPar.constAt(n, 4, 0)
0264                     << " inPar.constAt(n, 5, 0)=" << std::setprecision(9) << inPar.constAt(n, 5, 0)
0265                     << " inChg.constAt(n, 0, 0)=" << std::setprecision(9) << inChg.constAt(n, 0, 0));
0266     }
0267 #pragma omp simd
0268     for (int n = 0; n < NN; ++n) {
0269       dprint_np(n,
0270                 "propagation start, dump parameters"
0271                     << std::endl
0272                     << "pos = " << inPar.constAt(n, 0, 0) << " " << inPar.constAt(n, 1, 0) << " "
0273                     << inPar.constAt(n, 2, 0) << std::endl
0274                     << "mom (cart) = " << std::cos(inPar.constAt(n, 4, 0)) / inPar.constAt(n, 3, 0) << " "
0275                     << std::sin(inPar.constAt(n, 4, 0)) / inPar.constAt(n, 3, 0) << " "
0276                     << 1. / (inPar.constAt(n, 3, 0) * tan(inPar.constAt(n, 5, 0))) << " r="
0277                     << std::sqrt(inPar.constAt(n, 0, 0) * inPar.constAt(n, 0, 0) +
0278                                  inPar.constAt(n, 1, 0) * inPar.constAt(n, 1, 0))
0279                     << " pT=" << 1. / std::abs(inPar.constAt(n, 3, 0)) << " q=" << inChg.constAt(n, 0, 0)
0280                     << " targetZ=" << msZ.constAt(n, 0, 0) << std::endl);
0281     }
0282 
0283     float pt[NN];
0284 #pragma omp simd
0285     for (int n = 0; n < NN; ++n) {
0286       pt[n] = 1.f / ipt[n];
0287     }
0288 
0289     //no trig approx here, phi can be large
0290     float cosP[NN];
0291     float sinP[NN];
0292 #pragma omp simd
0293     for (int n = 0; n < NN; ++n) {
0294       cosP[n] = std::cos(phiin[n]);
0295     }
0296 
0297 #pragma omp simd
0298     for (int n = 0; n < NN; ++n) {
0299       sinP[n] = std::sin(phiin[n]);
0300     }
0301 
0302     float cosT[NN];
0303     float sinT[NN];
0304 #pragma omp simd
0305     for (int n = 0; n < NN; ++n) {
0306       cosT[n] = std::cos(theta[n]);
0307     }
0308 
0309 #pragma omp simd
0310     for (int n = 0; n < NN; ++n) {
0311       sinT[n] = std::sin(theta[n]);
0312     }
0313 
0314     float tanT[NN];
0315     float icos2T[NN];
0316     float pxin[NN];
0317     float pyin[NN];
0318 #pragma omp simd
0319     for (int n = 0; n < NN; ++n) {
0320       tanT[n] = sinT[n] / cosT[n];
0321       icos2T[n] = 1.f / (cosT[n] * cosT[n]);
0322       pxin[n] = cosP[n] * pt[n];
0323       pyin[n] = sinP[n] * pt[n];
0324     }
0325 
0326     float deltaZ[NN];
0327     float alpha[NN];
0328 #pragma omp simd
0329     for (int n = 0; n < NN; ++n) {
0330       deltaZ[n] = zout[n] - zin[n];
0331       alpha[n] = deltaZ[n] * tanT[n] * ipt[n] * kinv[n];
0332     }
0333 
0334     float cosahTmp[NN];
0335     float sinahTmp[NN];
0336     if constexpr (Config::useTrigApprox) {
0337 #if !defined(__INTEL_COMPILER)
0338 #pragma omp simd
0339 #endif
0340       for (int n = 0; n < NN; ++n) {
0341         sincos4(alpha[n] * 0.5f, sinahTmp[n], cosahTmp[n]);
0342       }
0343     } else {
0344 #if !defined(__INTEL_COMPILER)
0345 #pragma omp simd
0346 #endif
0347       for (int n = 0; n < NN; ++n) {
0348         cosahTmp[n] = std::cos(alpha[n] * 0.5f);
0349       }
0350 #if !defined(__INTEL_COMPILER)
0351 #pragma omp simd
0352 #endif
0353       for (int n = 0; n < NN; ++n) {
0354         sinahTmp[n] = std::sin(alpha[n] * 0.5f);
0355       }
0356     }
0357 
0358     float cosah[NN];
0359     float sinah[NN];
0360     float cosa[NN];
0361     float sina[NN];
0362 #pragma omp simd
0363     for (int n = 0; n < NN; ++n) {
0364       cosah[n] = cosahTmp[n];
0365       sinah[n] = sinahTmp[n];
0366       cosa[n] = 1.f - 2.f * sinah[n] * sinah[n];
0367       sina[n] = 2.f * sinah[n] * cosah[n];
0368     }
0369 
0370 //update parameters
0371 #pragma omp simd
0372     for (int n = 0; n < NN; ++n) {
0373       outPar.At(n, 0, 0) = outPar.At(n, 0, 0) + 2.f * k[n] * sinah[n] * (pxin[n] * cosah[n] - pyin[n] * sinah[n]);
0374       outPar.At(n, 1, 0) = outPar.At(n, 1, 0) + 2.f * k[n] * sinah[n] * (pyin[n] * cosah[n] + pxin[n] * sinah[n]);
0375       outPar.At(n, 2, 0) = zout[n];
0376       outPar.At(n, 4, 0) = phiin[n] + alpha[n];
0377     }
0378 
0379 #pragma omp simd
0380     for (int n = 0; n < NN; ++n) {
0381       dprint_np(n,
0382                 "propagation to Z end (OLD), dump parameters\n"
0383                     << "   pos = " << outPar(n, 0, 0) << " " << outPar(n, 1, 0) << " " << outPar(n, 2, 0) << "\t\t r="
0384                     << std::sqrt(outPar(n, 0, 0) * outPar(n, 0, 0) + outPar(n, 1, 0) * outPar(n, 1, 0)) << std::endl
0385                     << "   mom = " << outPar(n, 3, 0) << " " << outPar(n, 4, 0) << " " << outPar(n, 5, 0) << std::endl
0386                     << " cart= " << std::cos(outPar(n, 4, 0)) / outPar(n, 3, 0) << " "
0387                     << std::sin(outPar(n, 4, 0)) / outPar(n, 3, 0) << " "
0388                     << 1. / (outPar(n, 3, 0) * tan(outPar(n, 5, 0))) << "\t\tpT=" << 1. / std::abs(outPar(n, 3, 0))
0389                     << std::endl);
0390     }
0391 
0392     float pxcaMpysa[NN];
0393 #pragma omp simd
0394     for (int n = 0; n < NN; ++n) {
0395       pxcaMpysa[n] = pxin[n] * cosa[n] - pyin[n] * sina[n];
0396     }
0397 
0398 #pragma omp simd
0399     for (int n = 0; n < NN; ++n) {
0400       errorProp(n, 0, 2) = -tanT[n] * ipt[n] * pxcaMpysa[n];
0401       errorProp(n, 0, 3) =
0402           k[n] * pt[n] * pt[n] *
0403           (cosP[n] * (alpha[n] * cosa[n] - sina[n]) + sinP[n] * 2.f * sinah[n] * (sinah[n] - alpha[n] * cosah[n]));
0404       errorProp(n, 0, 4) = -2.f * k[n] * pt[n] * sinah[n] * (sinP[n] * cosah[n] + cosP[n] * sinah[n]);
0405       errorProp(n, 0, 5) = deltaZ[n] * ipt[n] * pxcaMpysa[n] * icos2T[n];
0406     }
0407 
0408     float pycaPpxsa[NN];
0409 #pragma omp simd
0410     for (int n = 0; n < NN; ++n) {
0411       pycaPpxsa[n] = pyin[n] * cosa[n] + pxin[n] * sina[n];
0412     }
0413 
0414 #pragma omp simd
0415     for (int n = 0; n < NN; ++n) {
0416       errorProp(n, 1, 2) = -tanT[n] * ipt[n] * pycaPpxsa[n];
0417       errorProp(n, 1, 3) =
0418           k[n] * pt[n] * pt[n] *
0419           (sinP[n] * (alpha[n] * cosa[n] - sina[n]) - cosP[n] * 2.f * sinah[n] * (sinah[n] - alpha[n] * cosah[n]));
0420       errorProp(n, 1, 4) = 2.f * k[n] * pt[n] * sinah[n] * (cosP[n] * cosah[n] - sinP[n] * sinah[n]);
0421       errorProp(n, 1, 5) = deltaZ[n] * ipt[n] * pycaPpxsa[n] * icos2T[n];
0422     }
0423 
0424 #pragma omp simd
0425     for (int n = 0; n < NN; ++n) {
0426       errorProp(n, 4, 2) = -ipt[n] * tanT[n] * kinv[n];
0427       errorProp(n, 4, 3) = tanT[n] * deltaZ[n] * kinv[n];
0428       errorProp(n, 4, 5) = ipt[n] * deltaZ[n] * kinv[n] * icos2T[n];
0429     }
0430 
0431 #pragma omp simd
0432     for (int n = 0; n < NN; ++n) {
0433       dprint_np(
0434           n,
0435           "propagation end, dump parameters"
0436               << std::endl
0437               << "pos = " << outPar.At(n, 0, 0) << " " << outPar.At(n, 1, 0) << " " << outPar.At(n, 2, 0) << std::endl
0438               << "mom (cart) = " << std::cos(outPar.At(n, 4, 0)) / outPar.At(n, 3, 0) << " "
0439               << std::sin(outPar.At(n, 4, 0)) / outPar.At(n, 3, 0) << " "
0440               << 1. / (outPar.At(n, 3, 0) * tan(outPar.At(n, 5, 0)))
0441               << " r=" << std::sqrt(outPar.At(n, 0, 0) * outPar.At(n, 0, 0) + outPar.At(n, 1, 0) * outPar.At(n, 1, 0))
0442               << " pT=" << 1. / std::abs(outPar.At(n, 3, 0)) << std::endl);
0443     }
0444 
0445     // PROP-FAIL-ENABLE Disabled to keep physics changes minimal.
0446     // To be reviewed, enabled and processed accordingly elsewhere.
0447     /*
0448     // Check for errors, set fail-flag.
0449     for (int n = 0; n < NN; ++n) {
0450       // We propagate for alpha: mark fail when prop angle more than pi/2
0451       if (std::abs(alpha[n]) > 1.57) {
0452         dprintf("helixAtZ: more than quarter turn, alpha = %f\n", alpha[n]);
0453         outFailFlag[n] = 1;
0454       } else {
0455         // Have we reached desired z? We can't know, we copy desired z to actual z.
0456         // Are we close to apex? Same condition as in propToR, 12.5 deg, cos(78.5deg) = 0.2
0457         float dotp = (outPar.At(n, 0, 0) * std::cos(outPar.At(n, 4, 0)) +
0458                       outPar.At(n, 1, 0) * std::sin(outPar.At(n, 4, 0))) /
0459                      hipo(outPar.At(n, 0, 0), outPar.At(n, 1, 0));
0460         if (dotp < 0.2 || dotp < 0) {
0461           dprintf("helixAtZ: dot product bad, dotp = %f\n", dotp);
0462           outFailFlag[n] = 1;
0463         }
0464       }
0465     }
0466     */
0467 
0468 #ifdef DEBUG
0469     if (debug && g_debug) {
0470       for (int n = 0; n < N_proc; ++n) {
0471         dmutex_guard;
0472         std::cout << n << ": jacobian" << std::endl;
0473         printf("%5f %5f %5f %5f %5f %5f\n",
0474                errorProp(n, 0, 0),
0475                errorProp(n, 0, 1),
0476                errorProp(n, 0, 2),
0477                errorProp(n, 0, 3),
0478                errorProp(n, 0, 4),
0479                errorProp(n, 0, 5));
0480         printf("%5f %5f %5f %5f %5f %5f\n",
0481                errorProp(n, 1, 0),
0482                errorProp(n, 1, 1),
0483                errorProp(n, 1, 2),
0484                errorProp(n, 1, 3),
0485                errorProp(n, 1, 4),
0486                errorProp(n, 1, 5));
0487         printf("%5f %5f %5f %5f %5f %5f\n",
0488                errorProp(n, 2, 0),
0489                errorProp(n, 2, 1),
0490                errorProp(n, 2, 2),
0491                errorProp(n, 2, 3),
0492                errorProp(n, 2, 4),
0493                errorProp(n, 2, 5));
0494         printf("%5f %5f %5f %5f %5f %5f\n",
0495                errorProp(n, 3, 0),
0496                errorProp(n, 3, 1),
0497                errorProp(n, 3, 2),
0498                errorProp(n, 3, 3),
0499                errorProp(n, 3, 4),
0500                errorProp(n, 3, 5));
0501         printf("%5f %5f %5f %5f %5f %5f\n",
0502                errorProp(n, 4, 0),
0503                errorProp(n, 4, 1),
0504                errorProp(n, 4, 2),
0505                errorProp(n, 4, 3),
0506                errorProp(n, 4, 4),
0507                errorProp(n, 4, 5));
0508         printf("%5f %5f %5f %5f %5f %5f\n",
0509                errorProp(n, 5, 0),
0510                errorProp(n, 5, 1),
0511                errorProp(n, 5, 2),
0512                errorProp(n, 5, 3),
0513                errorProp(n, 5, 4),
0514                errorProp(n, 5, 5));
0515       }
0516     }
0517 #endif
0518   }
0519 
0520 }  // namespace mkfit
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!