Remove the commented-out new implementation of prop-to-R that separates parameter and error propagation (as is done in prop-to-plane).

Author: osschar

RecoTracker/MkFitCore/src/PropagationMPlex.cc

@@ -488,19 +488,7 @@ namespace mkfit {
     errorProp.setVal(0.f);
     outFailFlag.setVal(0.f);
 
-    //helixAtRFromIterativeCCS_impl_new(inPar, inChg, msRad, outPar, errorProp, outFailFlag, 0, NN, N_proc, pflags);
     helixAtRFromIterativeCCS_impl(inPar, inChg, msRad, outPar, errorProp, outFailFlag, 0, NN, N_proc, pflags);
-    /*
-    //float nv = errorProp(0,0,0);
-
-    outPar = inPar;
-    errorProp.setVal(0.f);
-    outFailFlag.setVal(0.f);
-
-    helixAtRFromIterativeCCS_impl(inPar, inChg, msRad, outPar, errorProp, outFailFlag, 0, NN, N_proc, pflags);
-    //float ov = errorProp(0,0,0);
-    assert(0);
-    */
   }
 
   void propagateHelixToRMPlex(const MPlexLS& inErr,
@@ -1247,33 +1235,6 @@ namespace mkfit {
       }
     }
     // }
-
-    // This dump is now out of its place as similarity is done with matriplex ops.
-    /*
-#ifdef DEBUG
-   {
-     dmutex_guard;
-     for (int kk = 0; kk < N_proc; ++kk)
-     {
-       dprintf("outErr %d\n", kk);
-       for (int i = 0; i < 6; ++i) { for (int j = 0; j < 6; ++j)
-           dprintf("%8f ", outErr.At(kk,i,j)); printf("\n");
-       } dprintf("\n");
-
-       dprintf("outPar %d\n", kk);
-       for (int i = 0; i < 6; ++i) {
-           dprintf("%8f ", outPar.At(kk,i,0)); printf("\n");
-       } dprintf("\n");
-       if (std::abs(outPar.At(kk,2,0) - msZ.constAt(kk, 0, 0)) > 0.0001) {
-         float pt = 1.0f / inPar.constAt(kk,3,0);
-	 dprint_np(kk, "DID NOT GET TO Z, dZ=" << std::abs(outPar.At(kk,2,0) - msZ.constAt(kk, 0, 0))
-		   << " z=" << msZ.constAt(kk, 0, 0) << " zin=" << inPar.constAt(kk,2,0) << " zout=" << outPar.At(kk,2,0) << std::endl
-		   << "pt=" << pt << " pz=" << pt/std::tan(inPar.constAt(kk,5,0)));
-       }
-     }
-   }
-#endif
-   */
   }
 
   //==============================================================================

RecoTracker/MkFitCore/src/PropagationMPlex.icc

@@ -456,299 +456,6 @@ static inline void helixAtPlane_impl(const Tf& __restrict__ inPar,
   parsAndErrPropFromPathL_impl(inPar, inChg, outPar, kinv, s, errorProp, nmin, nmax, N_proc, pf);
 }
 
-/*
-// this function just calculates the path length (using the iterative approach as before)
-// and then calls parsAndErrPropFromPathL_impl for error propagation
-template <typename Tf, typename Ti, typename TfLL1, typename Tf11, typename TfLLL>
-static inline void helixAtRFromIterativeCCS_impl_new(const Tf& __restrict__ inPar,
-                                                     const Ti& __restrict__ inChg,
-                                                     const Tf11& __restrict__ msRad,
-                                                     TfLL1& __restrict__ outPar,
-                                                     TfLLL& __restrict__ errorProp,
-                                                     Ti& __restrict__ outFailFlag,  // expected to be initialized to 0
-                                                     const int nmin,
-                                                     const int nmax,
-                                                     const int N_proc,
-                                                     const PropagationFlags& pf) {
-
-#pragma omp simd
-  for (int n = nmin; n < nmax; ++n) {
-    //initialize erroProp to identity matrix
-    errorProp(n, 0, 0) = 1.f;
-    errorProp(n, 1, 1) = 1.f;
-    errorProp(n, 2, 2) = 1.f;
-    errorProp(n, 3, 3) = 1.f;
-    errorProp(n, 4, 4) = 1.f;
-    errorProp(n, 5, 5) = 1.f;
-  }
-  float r0[nmax - nmin];
-#pragma omp simd
-  for (int n = nmin; n < nmax; ++n) {
-    //initialize erroProp to identity matrix
-    r0[n - nmin] = hipo(inPar(n, 0, 0), inPar(n, 1, 0));
-  }
-  float k[nmax - nmin];
-  if (pf.use_param_b_field) {
-#pragma omp simd
-    for (int n = nmin; n < nmax; ++n) {
-      k[n - nmin] = inChg(n, 0, 0) * 100.f / (-Const::sol * Config::bFieldFromZR(inPar(n, 2, 0), r0[n - nmin]));
-    }
-  } else {
-#pragma omp simd
-    for (int n = nmin; n < nmax; ++n) {
-      k[n - nmin] = inChg(n, 0, 0) * 100.f / (-Const::sol * Config::Bfield);
-    }
-  }
-  float r[nmax - nmin];
-#pragma omp simd
-  for (int n = nmin; n < nmax; ++n) {
-    r[n - nmin] = msRad(n, 0, 0);
-  }
-  float xin[nmax - nmin];
-  float yin[nmax - nmin];
-  float ipt[nmax - nmin];
-  float phiin[nmax - nmin];
-#pragma omp simd
-  for (int n = nmin; n < nmax; ++n) {
-    xin[n - nmin] = inPar(n, 0, 0);
-    yin[n - nmin] = inPar(n, 1, 0);
-    ipt[n - nmin] = inPar(n, 3, 0);
-    phiin[n - nmin] = inPar(n, 4, 0);
-  }
-
-  for (int n = nmin; n < nmax; ++n) {
-    dprint_np(n,
-              "input parameters"
-                  << " inPar(n, 0, 0)=" << std::setprecision(9) << inPar(n, 0, 0) << " inPar(n, 1, 0)="
-                  << std::setprecision(9) << inPar(n, 1, 0) << " inPar(n, 2, 0)=" << std::setprecision(9)
-                  << inPar(n, 2, 0) << " inPar(n, 3, 0)=" << std::setprecision(9) << inPar(n, 3, 0)
-                  << " inPar(n, 4, 0)=" << std::setprecision(9) << inPar(n, 4, 0)
-                  << " inPar(n, 5, 0)=" << std::setprecision(9) << inPar(n, 5, 0));
-  }
-
-  float kinv[nmax - nmin];
-  float pt[nmax - nmin];
-#pragma omp simd
-  for (int n = nmin; n < nmax; ++n) {
-    kinv[n - nmin] = 1.f / k[n - nmin];
-    pt[n - nmin] = 1.f / ipt[n - nmin];
-  }
-  float D[nmax - nmin];
-  float cosa[nmax - nmin];
-  float sina[nmax - nmin];
-  float cosah[nmax - nmin];
-  float sinah[nmax - nmin];
-  float id[nmax - nmin];
-
-#pragma omp simd
-  for (int n = nmin; n < nmax; ++n) {
-    D[n - nmin] = 0.;
-  }
-
-  //no trig approx here, phi can be large
-  float cosPorT[nmax - nmin];
-  float sinPorT[nmax - nmin];
-#pragma omp simd
-  for (int n = nmin; n < nmax; ++n) {
-    cosPorT[n - nmin] = std::cos(phiin[n - nmin]);
-  }
-#pragma omp simd
-  for (int n = nmin; n < nmax; ++n) {
-    sinPorT[n - nmin] = std::sin(phiin[n - nmin]);
-  }
-
-  float pxin[nmax - nmin];
-  float pyin[nmax - nmin];
-#pragma omp simd
-  for (int n = nmin; n < nmax; ++n) {
-    pxin[n - nmin] = cosPorT[n - nmin] * pt[n - nmin];
-    pyin[n - nmin] = sinPorT[n - nmin] * pt[n - nmin];
-  }
-
-  for (int n = nmin; n < nmax; ++n) {
-    dprint_np(n,
-              "k=" << std::setprecision(9) << k[n - nmin] << " pxin=" << std::setprecision(9) << pxin[n - nmin]
-                   << " pyin=" << std::setprecision(9) << pyin[n - nmin] << " cosPorT=" << std::setprecision(9)
-                   << cosPorT[n - nmin] << " sinPorT=" << std::setprecision(9) << sinPorT[n - nmin]
-                   << " pt=" << std::setprecision(9) << pt[n - nmin]);
-  }
-
-  float oodotp[nmax - nmin];
-  float pxinold[nmax - nmin];
-
-  CMS_UNROLL_LOOP_COUNT(Config::Niter)
-  for (int i = 0; i < Config::Niter; ++i) {
-#pragma omp simd
-    for (int n = nmin; n < nmax; ++n) {
-      //compute distance and path for the current iteration
-      r0[n - nmin] = hipo(xin[n - nmin], yin[n - nmin]);
-    }
-
-    // Use one over dot product of transverse momentum and radial
-    // direction to scale the step. Propagation is prevented from reaching
-    // too close to the apex (dotp > 0.2).
-    // - Can / should we come up with a better approximation?
-    // - Can / should take +/- curvature into account?
-
-#pragma omp simd
-    for (int n = nmin; n < nmax; ++n) {
-      oodotp[n - nmin] =
-          r0[n - nmin] * pt[n - nmin] / (pxin[n - nmin] * xin[n - nmin] + pyin[n - nmin] * yin[n - nmin]);
-    }
-
-#pragma omp simd
-    for (int n = nmin; n < nmax; ++n) {
-      if (oodotp[n - nmin] > 5.0f || oodotp[n - nmin] < 0)  // 0.2 is 78.5 deg
-      {
-        outFailFlag(n, 0, 0) = 1;
-        oodotp[n - nmin] = 0.0f;
-      } else if (r[n - nmin] - r0[n - nmin] < 0.0f && pt[n - nmin] < 1.0f) {
-        // Scale down the correction for low-pT ingoing tracks.
-        oodotp[n - nmin] = 1.0f + (oodotp[n - nmin] - 1.0f) * pt[n - nmin];
-      }
-    }
-
-#pragma omp simd
-    for (int n = nmin; n < nmax; ++n) {
-      // Can we come up with a better approximation?
-      // Should take +/- curvature into account.
-      id[n - nmin] = (r[n - nmin] - r0[n - nmin]) * oodotp[n - nmin];
-    }
-
-#pragma omp simd
-    for (int n = nmin; n < nmax; ++n) {
-      D[n - nmin] += id[n - nmin];
-    }
-
-    if constexpr (Config::useTrigApprox) {
-#if !defined(__INTEL_COMPILER)
-#pragma omp simd
-#endif
-      for (int n = nmin; n < nmax; ++n) {
-        sincos4(id[n - nmin] * ipt[n - nmin] * kinv[n - nmin] * 0.5f, sinah[n - nmin], cosah[n - nmin]);
-      }
-    } else {
-#if !defined(__INTEL_COMPILER)
-#pragma omp simd
-#endif
-      for (int n = nmin; n < nmax; ++n) {
-        cosah[n - nmin] = std::cos(id[n - nmin] * ipt[n - nmin] * kinv[n - nmin] * 0.5f);
-        sinah[n - nmin] = std::sin(id[n - nmin] * ipt[n - nmin] * kinv[n - nmin] * 0.5f);
-      }
-    }
-
-#pragma omp simd
-    for (int n = nmin; n < nmax; ++n) {
-      cosa[n - nmin] = 1.f - 2.f * sinah[n - nmin] * sinah[n - nmin];
-      sina[n - nmin] = 2.f * sinah[n - nmin] * cosah[n - nmin];
-    }
-
-    for (int n = nmin; n < nmax; ++n) {
-      dprint_np(n,
-                "Attempt propagation from r="
-                    << r0[n - nmin] << " to r=" << r[n - nmin] << std::endl
-                    << "   x=" << xin[n - nmin] << " y=" << yin[n - nmin] << " z=" << inPar(n, 2, 0)
-                    << " px=" << pxin[n - nmin] << " py=" << pyin[n - nmin]
-                    << " pz=" << pt[n - nmin] / std::tan(inPar(n, 5, 0)) << " q=" << inChg(n, 0, 0) << std::endl
-                    << "   r=" << std::setprecision(9) << r[n - nmin] << " r0=" << std::setprecision(9) << r0[n - nmin]
-                    << " id=" << std::setprecision(9) << id[n - nmin] << " dr=" << std::setprecision(9)
-                    << r[n - nmin] - r0[n - nmin] << " cosa=" << cosa[n - nmin] << " sina=" << sina[n - nmin]
-                    << " dir_cos(rad,pT)=" << 1.0f / oodotp[n - nmin]);
-    }
-
-#pragma omp simd
-    for (int n = nmin; n < nmax; ++n) {
-      //update parameters
-      xin[n - nmin] = xin[n - nmin] + 2.f * k[n - nmin] * sinah[n - nmin] *
-                                              (pxin[n - nmin] * cosah[n - nmin] - pyin[n - nmin] * sinah[n - nmin]);
-      yin[n - nmin] = yin[n - nmin] + 2.f * k[n - nmin] * sinah[n - nmin] *
-                                              (pyin[n - nmin] * cosah[n - nmin] + pxin[n - nmin] * sinah[n - nmin]);
-      pxinold[n - nmin] = pxin[n - nmin];  //copy before overwriting
-      pxin[n - nmin] = pxin[n - nmin] * cosa[n - nmin] - pyin[n - nmin] * sina[n - nmin];
-      pyin[n - nmin] = pyin[n - nmin] * cosa[n - nmin] + pxinold[n - nmin] * sina[n - nmin];
-    }
-    for (int n = nmin; n < nmax; ++n) {
-      dprint_np(n,
-                "outPar(n, 0, 0)=" << outPar(n, 0, 0) << " outPar(n, 1, 0)=" << outPar(n, 1, 0)
-                                   << " pxin=" << pxin[n - nmin] << " pyin=" << pyin[n - nmin]);
-    }
-  }  // iteration loop
-
-  //float s[nmax - nmin];
-  MPlexQF s;
-  for (int n = nmin; n < nmax; ++n) {
-    //s[n - nmin] = D[n - nmin]/std::sin(inPar(n, 5, 0));
-    s(n, 0, 0) = D[n - nmin]/std::sin(inPar(n, 5, 0));
-  }
-  parsAndErrPropFromPathL_impl(inPar, inChg, outPar, kinv, s, errorProp, nmin, nmax, N_proc, pf);
-
-  for (int n = nmin; n < nmax; ++n) {
-    dprint_np(n,
-              "propagation to R end (NEW), dump parameters\n"
-                  << "   D = " << D[n - nmin] << " alpha = " << D[n - nmin] * inPar(n, 3, 0) * kinv[n - nmin] << " kinv = " << kinv[n - nmin] << std::endl
-                  << "   pos = " << outPar(n, 0, 0) << " " << outPar(n, 1, 0) << " " << outPar(n, 2, 0) << "\t\t r="
-                  << std::sqrt(outPar(n, 0, 0) * outPar(n, 0, 0) + outPar(n, 1, 0) * outPar(n, 1, 0)) << std::endl
-                  << "   mom = " << outPar(n, 3, 0) << " " << outPar(n, 4, 0) << " " << outPar(n, 5, 0) << std::endl
-		  << " cart= " << std::cos(outPar(n, 4, 0)) / outPar(n, 3, 0) << " "
-                  << std::sin(outPar(n, 4, 0)) / outPar(n, 3, 0) << " " << 1. / (outPar(n, 3, 0) * tan(outPar(n, 5, 0)))
-                  << "\t\tpT=" << 1. / std::abs(outPar(n, 3, 0)) << std::endl);
-  }
-
-#ifdef DEBUG
-  for (int n = nmin; n < nmax; ++n) {
-    if (debug && g_debug && n < N_proc) {
-      dmutex_guard;
-      std::cout << n << ": jacobian" << std::endl;
-      printf("%5f %5f %5f %5f %5f %5f\n",
-             errorProp(n, 0, 0),
-             errorProp(n, 0, 1),
-             errorProp(n, 0, 2),
-             errorProp(n, 0, 3),
-             errorProp(n, 0, 4),
-             errorProp(n, 0, 5));
-      printf("%5f %5f %5f %5f %5f %5f\n",
-             errorProp(n, 1, 0),
-             errorProp(n, 1, 1),
-             errorProp(n, 1, 2),
-             errorProp(n, 1, 3),
-             errorProp(n, 1, 4),
-             errorProp(n, 1, 5));
-      printf("%5f %5f %5f %5f %5f %5f\n",
-             errorProp(n, 2, 0),
-             errorProp(n, 2, 1),
-             errorProp(n, 2, 2),
-             errorProp(n, 2, 3),
-             errorProp(n, 2, 4),
-             errorProp(n, 2, 5));
-      printf("%5f %5f %5f %5f %5f %5f\n",
-             errorProp(n, 3, 0),
-             errorProp(n, 3, 1),
-             errorProp(n, 3, 2),
-             errorProp(n, 3, 3),
-             errorProp(n, 3, 4),
-             errorProp(n, 3, 5));
-      printf("%5f %5f %5f %5f %5f %5f\n",
-             errorProp(n, 4, 0),
-             errorProp(n, 4, 1),
-             errorProp(n, 4, 2),
-             errorProp(n, 4, 3),
-             errorProp(n, 4, 4),
-             errorProp(n, 4, 5));
-      printf("%5f %5f %5f %5f %5f %5f\n",
-             errorProp(n, 5, 0),
-             errorProp(n, 5, 1),
-             errorProp(n, 5, 2),
-             errorProp(n, 5, 3),
-             errorProp(n, 5, 4),
-             errorProp(n, 5, 5));
-      printf("\n");
-    }
-  }
-#endif
-
-}
-*/
-
 template <typename Tf, typename Ti, typename TfLL1, typename Tf11, typename TfLLL>
 static inline void helixAtRFromIterativeCCS_impl(const Tf& __restrict__ inPar,
                                                  const Ti& __restrict__ inChg,