File indexing completed on 2024-10-17 22:59:07
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
0009 #include "Debug.h"
0010
0011 namespace {
0012 using namespace mkfit;
0013
0014 void MultHelixPropEndcap(const MPlexLL& A, const MPlexLS& B, MPlexLL& C) {
0015
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
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 }
0047
0048
0049
0050 namespace {}
0051
0052
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
0068
0069 outErr = inErr;
0070 outPar = inPar;
0071
0072 MPlexLL errorProp;
0073
0074
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
0120 MPlexLL temp;
0121 MultHelixPropEndcap(errorProp, outErr, temp);
0122 MultHelixPropTranspEndcap(errorProp, temp, outErr);
0123
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);
0183
0184
0185
0186
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
0208 #pragma omp simd
0209 for (int n = 0; n < NN; ++n) {
0210
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
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
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
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
0444
0445
0446
0447
0448
0449
0450
0451
0452
0453
0454
0455
0456
0457
0458
0459
0460
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 }