diff --git a/RecoTracker/MkFitCore/src/PropagationMPlex.cc b/RecoTracker/MkFitCore/src/PropagationMPlex.cc
index bc9e5328232a0..f6a14485a7acf 100644
--- a/RecoTracker/MkFitCore/src/PropagationMPlex.cc
+++ b/RecoTracker/MkFitCore/src/PropagationMPlex.cc
@@ -588,16 +588,16 @@ namespace mkfit {
   }
 
   //==============================================================================
-  
-  void __attribute__ ((optimize("no-inline"))) propagateHelixToZMPlex(const MPlexLS& inErr,
-                              const MPlexLV& inPar,
-                              const MPlexQI& inChg,
-                              const MPlexQF& msZ,
-                              MPlexLS& outErr,
-                              MPlexLV& outPar,
-                              const int N_proc,
-                              const PropagationFlags pflags,
-                              const MPlexQI* noMatEffPtr) {
+
+  void __attribute__((optimize("no-inline"))) propagateHelixToZMPlex(const MPlexLS& inErr,
+                                                                     const MPlexLV& inPar,
+                                                                     const MPlexQI& inChg,
+                                                                     const MPlexQF& msZ,
+                                                                     MPlexLS& outErr,
+                                                                     MPlexLV& outPar,
+                                                                     const int N_proc,
+                                                                     const PropagationFlags pflags,
+                                                                     const MPlexQI* noMatEffPtr) {
     // debug = true;
 
     outErr = inErr;
@@ -708,7 +708,7 @@ namespace mkfit {
       errorProp(n, 3, 3) = 1.f;
       errorProp(n, 4, 4) = 1.f;
       errorProp(n, 5, 5) = 1.f;
-}
+    }
     float zout[NN];
     float zin[NN];
     float ipt[NN];
@@ -728,14 +728,13 @@ namespace mkfit {
     if (pflags.use_param_b_field) {
 #pragma omp simd
       for (int n = 0; n < NN; ++n) {
-          k[n] = inChg.constAt(n, 0, 0) * 100.f / (-Const::sol * 
-                             Config::bFieldFromZR(zin[n], hipo(inPar.constAt(n, 0, 0), inPar.constAt(n, 1, 0))));
+        k[n] = inChg.constAt(n, 0, 0) * 100.f /
+               (-Const::sol * Config::bFieldFromZR(zin[n], hipo(inPar.constAt(n, 0, 0), inPar.constAt(n, 1, 0))));
       }
-    } else { 
+    } else {
 #pragma omp simd
       for (int n = 0; n < NN; ++n) {
-        k[n] = inChg.constAt(n, 0, 0) * 100.f / (-Const::sol * 
-                            Config::Bfield);
+        k[n] = inChg.constAt(n, 0, 0) * 100.f / (-Const::sol * Config::Bfield);
       }
     }
 
@@ -763,8 +762,8 @@ namespace mkfit {
     for (int n = 0; n < NN; ++n) {
       pt[n] = 1.f / ipt[n];
     }
-    
-      //no trig approx here, phi can be large
+
+    //no trig approx here, phi can be large
     float cosP[NN];
     float sinP[NN];
 #pragma omp simd
@@ -802,7 +801,6 @@ namespace mkfit {
     }
 #pragma omp simd
     for (int n = 0; n < NN; ++n) {
-
       //fixme, make this printout useful for propagation to z
       dprint_np(n,
                 std::endl
@@ -847,7 +845,7 @@ namespace mkfit {
       cosa[n] = 1.f - 2.f * sinah[n] * sinah[n];
       sina[n] = 2.f * sinah[n] * cosah[n];
     }
-      //update parameters
+    //update parameters
 #pragma omp simd
     for (int n = 0; n < NN; ++n) {
       outPar.At(n, 0, 0) = outPar.At(n, 0, 0) + 2.f * k[n] * sinah[n] * (pxin[n] * cosah[n] - pyin[n] * sinah[n]);
@@ -873,7 +871,9 @@ namespace mkfit {
 #pragma omp simd
     for (int n = 0; n < NN; ++n) {
       errorProp(n, 0, 2) = -tanT[n] * ipt[n] * pxcaMpysa[n];
-      errorProp(n, 0, 3) = k[n] * pt[n] * pt[n] * (cosP[n] * (alpha[n] * cosa[n] - sina[n]) + sinP[n] * 2.f * sinah[n] * (sinah[n] - alpha[n] * cosah[n]));
+      errorProp(n, 0, 3) =
+          k[n] * pt[n] * pt[n] *
+          (cosP[n] * (alpha[n] * cosa[n] - sina[n]) + sinP[n] * 2.f * sinah[n] * (sinah[n] - alpha[n] * cosah[n]));
       errorProp(n, 0, 4) = -2.f * k[n] * pt[n] * sinah[n] * (sinP[n] * cosah[n] + cosP[n] * sinah[n]);
       errorProp(n, 0, 5) = deltaZ[n] * ipt[n] * pxcaMpysa[n] * icos2T[n];
     }
@@ -887,7 +887,9 @@ namespace mkfit {
 #pragma omp simd
     for (int n = 0; n < NN; ++n) {
       errorProp(n, 1, 2) = -tanT[n] * ipt[n] * pycaPpxsa[n];
-      errorProp(n, 1, 3) = k[n] * pt[n] * pt[n] * (sinP[n] * (alpha[n] * cosa[n] - sina[n]) - cosP[n] * 2.f * sinah[n] * (sinah[n] - alpha[n] * cosah[n]));
+      errorProp(n, 1, 3) =
+          k[n] * pt[n] * pt[n] *
+          (sinP[n] * (alpha[n] * cosa[n] - sina[n]) - cosP[n] * 2.f * sinah[n] * (sinah[n] - alpha[n] * cosah[n]));
       errorProp(n, 1, 4) = 2.f * k[n] * pt[n] * sinah[n] * (cosP[n] * cosah[n] - sinP[n] * sinah[n]);
       errorProp(n, 1, 5) = deltaZ[n] * ipt[n] * pycaPpxsa[n] * icos2T[n];
     }
@@ -975,7 +977,6 @@ namespace mkfit {
                             MPlexLV& outPar,
                             const int N_proc,
                             const bool isBarrel) {
-
     float radL[NN];
 #pragma omp simd
     for (int n = 0; n < NN; ++n) {
@@ -984,25 +985,29 @@ namespace mkfit {
     float theta[NN];
 #pragma omp simd
     for (int n = 0; n < NN; ++n) {
-      if (radL[n] < 1e-13f) continue;  //ugly, please fixme
+      if (radL[n] < 1e-13f)
+        continue;  //ugly, please fixme
       theta[n] = outPar.constAt(n, 5, 0);
     }
     float pt[NN];
 #pragma omp simd
     for (int n = 0; n < NN; ++n) {
-      if (radL[n] < 1e-13f) continue;  //ugly, please fixme
+      if (radL[n] < 1e-13f)
+        continue;                             //ugly, please fixme
       pt[n] = 1.f / outPar.constAt(n, 3, 0);  //fixme, make sure it is positive?
     }
     float p[NN];
 #pragma omp simd
     for (int n = 0; n < NN; ++n) {
-      if (radL[n] < 1e-13f) continue;  //ugly, please fixme
+      if (radL[n] < 1e-13f)
+        continue;  //ugly, please fixme
       p[n] = pt[n] / std::sin(theta[n]);
     }
     float p2[NN];
 #pragma omp simd
     for (int n = 0; n < NN; ++n) {
-      if (radL[n] < 1e-13f) continue;  //ugly, please fixme
+      if (radL[n] < 1e-13f)
+        continue;  //ugly, please fixme
       p2[n] = p[n] * p[n];
     }
     constexpr float mpi = 0.140;       // m=140 MeV, pion
@@ -1010,117 +1015,133 @@ namespace mkfit {
     float beta2[NN];
 #pragma omp simd
     for (int n = 0; n < NN; ++n) {
-      if (radL[n] < 1e-13f) continue;  //ugly, please fixme
+      if (radL[n] < 1e-13f)
+        continue;  //ugly, please fixme
       beta2[n] = p2[n] / (p2[n] + mpi2);
     }
     float beta[NN];
 #pragma omp simd
     for (int n = 0; n < NN; ++n) {
-      if (radL[n] < 1e-13f) continue;  //ugly, please fixme
+      if (radL[n] < 1e-13f)
+        continue;  //ugly, please fixme
       beta[n] = std::sqrt(beta2[n]);
     }
-      //radiation lenght, corrected for the crossing angle (cos alpha from dot product of radius vector and momentum)
+    //radiation lenght, corrected for the crossing angle (cos alpha from dot product of radius vector and momentum)
     float invCos[NN];
 #pragma omp simd
     for (int n = 0; n < NN; ++n) {
-      if (radL[n] < 1e-13f) continue;  //ugly, please fixme
+      if (radL[n] < 1e-13f)
+        continue;  //ugly, please fixme
       invCos[n] = (isBarrel ? p[n] / pt[n] : 1.f / std::abs(std::cos(theta[n])));
     }
 #pragma omp simd
     for (int n = 0; n < NN; ++n) {
-      if (radL[n] < 1e-13f) continue;  //ugly, please fixme
+      if (radL[n] < 1e-13f)
+        continue;                     //ugly, please fixme
       radL[n] = radL[n] * invCos[n];  //fixme works only for barrel geom
     }
-      // multiple scattering
-      //vary independently phi and theta by the rms of the planar multiple scattering angle
-      // XXX-KMD radL < 0, see your fixme above! Repeating bailout
-//      if (radL < 1e-13f)
-//        continue;
-      // const float thetaMSC = 0.0136f*std::sqrt(radL)*(1.f+0.038f*std::log(radL))/(beta*p);// eq 32.15
-      // const float thetaMSC2 = thetaMSC*thetaMSC;
+    // multiple scattering
+    //vary independently phi and theta by the rms of the planar multiple scattering angle
+    // XXX-KMD radL < 0, see your fixme above! Repeating bailout
+    //      if (radL < 1e-13f)
+    //        continue;
+    // const float thetaMSC = 0.0136f*std::sqrt(radL)*(1.f+0.038f*std::log(radL))/(beta*p);// eq 32.15
+    // const float thetaMSC2 = thetaMSC*thetaMSC;
     float thetaMSC[NN];
 #pragma omp simd
     for (int n = 0; n < NN; ++n) {
-      if (radL[n] < 1e-13f) continue;  //ugly, please fixme
+      if (radL[n] < 1e-13f)
+        continue;                                                                     //ugly, please fixme
       thetaMSC[n] = 0.0136f * (1.f + 0.038f * std::log(radL[n])) / (beta[n] * p[n]);  // eq 32.15
     }
     float thetaMSC2[NN];
 #pragma omp simd
     for (int n = 0; n < NN; ++n) {
-      if (radL[n] < 1e-13f) continue;  //ugly, please fixme
+      if (radL[n] < 1e-13f)
+        continue;  //ugly, please fixme
       thetaMSC2[n] = thetaMSC[n] * thetaMSC[n] * radL[n];
     }
 #pragma omp simd
     for (int n = 0; n < NN; ++n) {
-      if (radL[n] < 1e-13f) continue;  //ugly, please fixme
+      if (radL[n] < 1e-13f)
+        continue;  //ugly, please fixme
       outErr.At(n, 4, 4) += thetaMSC2[n];
     }
-      // outErr.At(n, 4, 5) += thetaMSC2;
+    // outErr.At(n, 4, 5) += thetaMSC2;
 #pragma omp simd
     for (int n = 0; n < NN; ++n) {
-      if (radL[n] < 1e-13f) continue;  //ugly, please fixme
+      if (radL[n] < 1e-13f)
+        continue;  //ugly, please fixme
       outErr.At(n, 5, 5) += thetaMSC2[n];
     }
-      //std::cout << "beta=" << beta << " p=" << p << std::endl;
-      //std::cout << "multiple scattering thetaMSC=" << thetaMSC << " thetaMSC2=" << thetaMSC2 << " radL=" << radL << std::endl;
-      // energy loss
-      // XXX-KMD beta2 = 1 => 1 / sqrt(0)
-      // const float gamma = 1.f/std::sqrt(1.f - std::min(beta2, 0.999999f));
-      // const float gamma2 = gamma*gamma;
+    //std::cout << "beta=" << beta << " p=" << p << std::endl;
+    //std::cout << "multiple scattering thetaMSC=" << thetaMSC << " thetaMSC2=" << thetaMSC2 << " radL=" << radL << std::endl;
+    // energy loss
+    // XXX-KMD beta2 = 1 => 1 / sqrt(0)
+    // const float gamma = 1.f/std::sqrt(1.f - std::min(beta2, 0.999999f));
+    // const float gamma2 = gamma*gamma;
     float gamma2[NN];
 #pragma omp simd
     for (int n = 0; n < NN; ++n) {
-      if (radL[n] < 1e-13f) continue;  //ugly, please fixme
+      if (radL[n] < 1e-13f)
+        continue;  //ugly, please fixme
       gamma2[n] = (p2[n] + mpi2) / mpi2;
     }
     float gamma[NN];
 #pragma omp simd
     for (int n = 0; n < NN; ++n) {
-      if (radL[n] < 1e-13f) continue;  //ugly, please fixme
-      gamma[n]= std::sqrt(gamma2[n]);  //1.f/std::sqrt(1.f - std::min(beta2, 0.999999f));
+      if (radL[n] < 1e-13f)
+        continue;                       //ugly, please fixme
+      gamma[n] = std::sqrt(gamma2[n]);  //1.f/std::sqrt(1.f - std::min(beta2, 0.999999f));
     }
-    constexpr float me = 0.0005;            // m=0.5 MeV, electron
+    constexpr float me = 0.0005;  // m=0.5 MeV, electron
     float wmax[NN];
 #pragma omp simd
     for (int n = 0; n < NN; ++n) {
-      if (radL[n] < 1e-13f) continue;  //ugly, please fixme
+      if (radL[n] < 1e-13f)
+        continue;  //ugly, please fixme
       wmax[n] = 2.f * me * beta2[n] * gamma2[n] / (1.f + 2.f * gamma[n] * me / mpi + me * me / (mpi * mpi));
     }
     constexpr float I = 16.0e-9 * 10.75;
     float deltahalf[NN];
 #pragma omp simd
     for (int n = 0; n < NN; ++n) {
-      if (radL[n] < 1e-13f) continue;  //ugly, please fixme
+      if (radL[n] < 1e-13f)
+        continue;  //ugly, please fixme
       deltahalf[n] = std::log(28.816e-9f * std::sqrt(2.33f * 0.498f) / I) + std::log(beta[n] * gamma[n]) - 0.5f;
     }
     float dEdx[NN];
     float dedxtmp[NN];
 #pragma omp simd
     for (int n = 0; n < NN; ++n) {
-      if (radL[n] < 1e-13f) continue;  //ugly, please fixme
-      dedxtmp[n] = 2.f * (hitsXi.constAt(n, 0, 0) * invCos[n] *
-                        (0.5f * std::log(2.f * me * beta2[n] * gamma2[n] * wmax[n] / (I * I)) - beta2[n] - deltahalf[n]) / beta2[n]);
-      dEdx[n] =
-          beta[n] < 1.f ? dedxtmp[n]
-              : 0.f;  //protect against infs and nans
+      if (radL[n] < 1e-13f)
+        continue;  //ugly, please fixme
+      dedxtmp[n] =
+          2.f *
+          (hitsXi.constAt(n, 0, 0) * invCos[n] *
+           (0.5f * std::log(2.f * me * beta2[n] * gamma2[n] * wmax[n] / (I * I)) - beta2[n] - deltahalf[n]) / beta2[n]);
+      dEdx[n] = beta[n] < 1.f ? dedxtmp[n] : 0.f;  //protect against infs and nans
     }
-      // dEdx = dEdx*2.;//xi in cmssw is defined with an extra factor 0.5 with respect to formula 27.1 in pdg
-      //std::cout << "dEdx=" << dEdx << " delta=" << deltahalf << " wmax=" << wmax << " Xi=" << hitsXi.constAt(n,0,0) << std::endl;
+    // dEdx = dEdx*2.;//xi in cmssw is defined with an extra factor 0.5 with respect to formula 27.1 in pdg
+    //std::cout << "dEdx=" << dEdx << " delta=" << deltahalf << " wmax=" << wmax << " Xi=" << hitsXi.constAt(n,0,0) << std::endl;
     float dP[NN];
 #pragma omp simd
     for (int n = 0; n < NN; ++n) {
-      if (radL[n] < 1e-13f) continue;  //ugly, please fixme
-      dP[n]= propSign.constAt(n, 0, 0) * dEdx[n] / beta[n];
+      if (radL[n] < 1e-13f)
+        continue;  //ugly, please fixme
+      dP[n] = propSign.constAt(n, 0, 0) * dEdx[n] / beta[n];
     }
 #pragma omp simd
     for (int n = 0; n < NN; ++n) {
-      if (radL[n] < 1e-13f) continue;  //ugly, please fixme
+      if (radL[n] < 1e-13f)
+        continue;                                                            //ugly, please fixme
       outPar.At(n, 3, 0) = p[n] / (std::max(p[n] + dP[n], 0.001f) * pt[n]);  //stay above 1MeV
       //assume 100% uncertainty
     }
 #pragma omp simd
     for (int n = 0; n < NN; ++n) {
-      if (radL[n] < 1e-13f) continue;  //ugly, please fixme
+      if (radL[n] < 1e-13f)
+        continue;  //ugly, please fixme
       outErr.At(n, 3, 3) += dP[n] * dP[n] / (p2[n] * pt[n] * pt[n]);
     }
   }
diff --git a/RecoTracker/MkFitCore/src/PropagationMPlex.icc b/RecoTracker/MkFitCore/src/PropagationMPlex.icc
index 262d2de917b83..3c9d66db1ae73 100644
--- a/RecoTracker/MkFitCore/src/PropagationMPlex.icc
+++ b/RecoTracker/MkFitCore/src/PropagationMPlex.icc
@@ -25,36 +25,34 @@ static inline void helixAtRFromIterativeCCS_impl(const Tf& __restrict__ inPar,
     errorProp(n, 4, 4) = 1.f;
     errorProp(n, 5, 5) = 1.f;
   }
-  float r0[nmax-nmin];
+  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));
+    r0[n - nmin] = hipo(inPar(n, 0, 0), inPar(n, 1, 0));
   }
-  float k[nmax-nmin];
+  float k[nmax - nmin];
   if (pf.use_param_b_field) {
-    #pragma omp simd
+#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]));
+      k[n - nmin] = inChg(n, 0, 0) * 100.f / (-Const::sol * Config::bFieldFromZR(inPar(n, 2, 0), r0[n - nmin]));
     }
   } else {
-    #pragma omp simd
+#pragma omp simd
     for (int n = nmin; n < nmax; ++n) {
-      k[n-nmin] = inChg(n, 0, 0) * 100.f /
-                      (-Const::sol * Config::Bfield);
+      k[n - nmin] = inChg(n, 0, 0) * 100.f / (-Const::sol * Config::Bfield);
     }
   }
-  float r[nmax-nmin];
+  float r[nmax - nmin];
 #pragma omp simd
   for (int n = nmin; n < nmax; ++n) {
-    r[n-nmin] = msRad(n, 0, 0);
+    r[n - nmin] = msRad(n, 0, 0);
   }
-  float xin[nmax-nmin];
-  float yin[nmax-nmin];
-  float ipt[nmax-nmin];
-  float phiin[nmax-nmin];
-  float theta[nmax-nmin];
+  float xin[nmax - nmin];
+  float yin[nmax - nmin];
+  float ipt[nmax - nmin];
+  float phiin[nmax - nmin];
+  float theta[nmax - nmin];
 #pragma omp simd
   for (int n = nmin; n < nmax; ++n) {
     // if (std::abs(r-r0)<0.0001f) {
@@ -62,11 +60,11 @@ static inline void helixAtRFromIterativeCCS_impl(const Tf& __restrict__ inPar,
     // 	continue;
     // }
 
-    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);
-    theta[n-nmin] = inPar(n, 5, 0);
+    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);
+    theta[n - nmin] = inPar(n, 5, 0);
 
     dprint(std::endl);
   }
@@ -81,298 +79,332 @@ static inline void helixAtRFromIterativeCCS_impl(const Tf& __restrict__ inPar,
                   << " inPar(n, 5, 0)=" << std::setprecision(9) << inPar(n, 5, 0));
   }
 
-  float kinv[nmax-nmin];
-  float pt[nmax-nmin];
+  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];
+    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];
-
-    //no trig approx here, phi can be large
-  float cosPorT[nmax-nmin];
-  float sinPorT[nmax-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];
+
+  //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]);
+    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]);
+    sinPorT[n - nmin] = std::sin(phiin[n - nmin]);
   }
 
-  float pxin[nmax-nmin];
-  float pyin[nmax-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];
+    pxin[n - nmin] = cosPorT[n - nmin] * pt[n - nmin];
+    pyin[n - nmin] = sinPorT[n - nmin] * pt[n - nmin];
   }
 #pragma omp simd
   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]);
+              "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 dDdx[nmax-nmin];
-  float dDdy[nmax-nmin];
-  float dDdipt[nmax-nmin];
-  float dDdphi[nmax-nmin];
+  float dDdx[nmax - nmin];
+  float dDdy[nmax - nmin];
+  float dDdipt[nmax - nmin];
+  float dDdphi[nmax - nmin];
 
 #pragma omp simd
   for (int n = nmin; n < nmax; ++n) {
     //derivatives initialized to value for first iteration, i.e. distance = r-r0in
-    dDdx[n-nmin] = r0[n-nmin] > 0.f ? -xin[n-nmin] / r0[n-nmin] : 0.f;
+    dDdx[n - nmin] = r0[n - nmin] > 0.f ? -xin[n - nmin] / r0[n - nmin] : 0.f;
   }
 
 #pragma omp simd
   for (int n = nmin; n < nmax; ++n) {
-    dDdy[n-nmin] = r0[n-nmin] > 0.f ? -yin[n-nmin] / r0[n-nmin] : 0.f;
+    dDdy[n - nmin] = r0[n - nmin] > 0.f ? -yin[n - nmin] / r0[n - nmin] : 0.f;
   }
 
-  float oodotp[nmax-nmin];
-  float x[nmax-nmin];
-  float y[nmax-nmin];
-  float oor0[nmax-nmin];
-  float dadipt[nmax-nmin];
-  float dadx[nmax-nmin];
-  float dady[nmax-nmin];
-  float pxca[nmax-nmin];
-  float pxsa[nmax-nmin];
-  float pyca[nmax-nmin];
-  float pysa[nmax-nmin];
-  float tmp[nmax-nmin];
-  float tmpx[nmax-nmin];
-  float tmpy[nmax-nmin];
-  float pxinold[nmax-nmin]; 
+  float oodotp[nmax - nmin];
+  float x[nmax - nmin];
+  float y[nmax - nmin];
+  float oor0[nmax - nmin];
+  float dadipt[nmax - nmin];
+  float dadx[nmax - nmin];
+  float dady[nmax - nmin];
+  float pxca[nmax - nmin];
+  float pxsa[nmax - nmin];
+  float pyca[nmax - nmin];
+  float pysa[nmax - nmin];
+  float tmp[nmax - nmin];
+  float tmpx[nmax - nmin];
+  float tmpy[nmax - nmin];
+  float pxinold[nmax - nmin];
 
   for (int i = 0; i < Config::Niter; ++i) {
-
-    #pragma omp simd
+#pragma omp simd
     for (int n = nmin; n < nmax; ++n) {
       //compute distance and path for the current iteration
-      r0[n-nmin] = hipo(outPar(n, 0, 0), outPar(n, 1, 0));
+      r0[n - nmin] = hipo(outPar(n, 0, 0), outPar(n, 1, 0));
     }
 
-      // Use one over dot produce 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?
+    // Use one over dot produce 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
+#pragma omp simd
     for (int n = nmin; n < nmax; ++n) {
-      oodotp[n-nmin] = r0[n-nmin] * pt[n-nmin] / (pxin[n-nmin] * outPar(n, 0, 0) + pyin[n-nmin] * outPar(n, 1, 0));
+      oodotp[n - nmin] =
+          r0[n - nmin] * pt[n - nmin] / (pxin[n - nmin] * outPar(n, 0, 0) + pyin[n - nmin] * outPar(n, 1, 0));
     }
 
-    #pragma omp simd
+#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
+      if (oodotp[n - nmin] > 5.0f || oodotp[n - nmin] < 0)  // 0.2 is 78.5 deg
       {
         outFailFlag(n, 0, 0) = 1;
       }
-    } 
+    }
 
-    #pragma omp simd
+#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] = (oodotp[n-nmin] > 5.0f || oodotp[n-nmin] < 0) ? 0.0f : (r[n-nmin] - r0[n-nmin]) * oodotp[n-nmin];
+      // Can we come up with a better approximation?
+      // Should take +/- curvature into account.
+      id[n - nmin] =
+          (oodotp[n - nmin] > 5.0f || oodotp[n - nmin] < 0) ? 0.0f : (r[n - nmin] - r0[n - nmin]) * oodotp[n - nmin];
     }
- 
-    #pragma omp simd
+
+#pragma omp simd
     for (int n = nmin; n < nmax; ++n) {
-      D[n-nmin] += id[n-nmin];
+      D[n - nmin] += id[n - nmin];
     }
 
     if constexpr (Config::useTrigApprox) {
-      #pragma omp simd
+#pragma omp simd
       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]);
+        sincos4(id[n - nmin] * ipt[n - nmin] * kinv[n - nmin] * 0.5f, sinah[n - nmin], cosah[n - nmin]);
       }
     } else {
-      #pragma omp simd
+#pragma omp simd
       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);
+        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
+#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];
+      cosa[n - nmin] = 1.f - 2.f * sinah[n - nmin] * sinah[n - nmin];
+      sina[n - nmin] = 2.f * sinah[n - nmin] * cosah[n - nmin];
     }
 
-    #pragma omp simd
+#pragma omp simd
     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(theta[n-nmin]) << " 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]);
+                    << 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(theta[n - nmin]) << " 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]);
     }
 
     //update derivatives on total distance
     if (i + 1 != Config::Niter) {
-      #pragma omp simd
+#pragma omp simd
       for (int n = nmin; n < nmax; ++n) {
-        x[n-nmin] = outPar(n, 0, 0);
-        y[n-nmin] = outPar(n, 1, 0);
+        x[n - nmin] = outPar(n, 0, 0);
+        y[n - nmin] = outPar(n, 1, 0);
       }
-      #pragma omp simd
+#pragma omp simd
       for (int n = nmin; n < nmax; ++n) {
-        oor0[n-nmin]  = (r0[n-nmin] > 0.f && std::abs(r[n-nmin] - r0[n-nmin]) < 0.0001f) ? 1.f / r0[n-nmin] : 0.f;
+        oor0[n - nmin] =
+            (r0[n - nmin] > 0.f && std::abs(r[n - nmin] - r0[n - nmin]) < 0.0001f) ? 1.f / r0[n - nmin] : 0.f;
       }
-      #pragma omp simd
+#pragma omp simd
       for (int n = nmin; n < nmax; ++n) {
-        dadipt[n-nmin]  = id[n-nmin] * kinv[n-nmin];
-        dadx[n-nmin]  = -x[n-nmin] * ipt[n-nmin] * kinv[n-nmin] * oor0[n-nmin];
-        dady[n-nmin]  = -y[n-nmin] * ipt[n-nmin] * kinv[n-nmin] * oor0[n-nmin];
-        pxca[n-nmin] = pxin[n-nmin] * cosa[n-nmin];
-        pxsa[n-nmin] = pxin[n-nmin] * sina[n-nmin];
-        pyca[n-nmin] = pyin[n-nmin] * cosa[n-nmin];
-        pysa[n-nmin] = pyin[n-nmin] * sina[n-nmin];
-        tmpx[n-nmin] = k[n-nmin] * dadx[n-nmin];
+        dadipt[n - nmin] = id[n - nmin] * kinv[n - nmin];
+        dadx[n - nmin] = -x[n - nmin] * ipt[n - nmin] * kinv[n - nmin] * oor0[n - nmin];
+        dady[n - nmin] = -y[n - nmin] * ipt[n - nmin] * kinv[n - nmin] * oor0[n - nmin];
+        pxca[n - nmin] = pxin[n - nmin] * cosa[n - nmin];
+        pxsa[n - nmin] = pxin[n - nmin] * sina[n - nmin];
+        pyca[n - nmin] = pyin[n - nmin] * cosa[n - nmin];
+        pysa[n - nmin] = pyin[n - nmin] * sina[n - nmin];
+        tmpx[n - nmin] = k[n - nmin] * dadx[n - nmin];
       }
 
-      #pragma omp simd
+#pragma omp simd
       for (int n = nmin; n < nmax; ++n) {
-        dDdx[n-nmin] -= (x[n-nmin] * (1.f + tmpx[n-nmin] * (pxca[n-nmin] - pysa[n-nmin])) + y[n-nmin] * tmpx[n-nmin] * (pyca[n-nmin] + pxsa[n-nmin])) * oor0[n-nmin];
+        dDdx[n - nmin] -= (x[n - nmin] * (1.f + tmpx[n - nmin] * (pxca[n - nmin] - pysa[n - nmin])) +
+                           y[n - nmin] * tmpx[n - nmin] * (pyca[n - nmin] + pxsa[n - nmin])) *
+                          oor0[n - nmin];
       }
 
-      #pragma omp simd
+#pragma omp simd
       for (int n = nmin; n < nmax; ++n) {
-        tmpy[n-nmin] = k[n-nmin] * dady[n-nmin];
+        tmpy[n - nmin] = k[n - nmin] * dady[n - nmin];
       }
-      #pragma omp simd
+#pragma omp simd
       for (int n = nmin; n < nmax; ++n) {
-        dDdy[n-nmin] -= (x[n-nmin] * tmpy[n-nmin] * (pxca[n-nmin] - pysa[n-nmin]) + y[n-nmin] * (1.f + tmpy[n-nmin] * (pyca[n-nmin] + pxsa[n-nmin]))) * oor0[n-nmin];
+        dDdy[n - nmin] -= (x[n - nmin] * tmpy[n - nmin] * (pxca[n - nmin] - pysa[n - nmin]) +
+                           y[n - nmin] * (1.f + tmpy[n - nmin] * (pyca[n - nmin] + pxsa[n - nmin]))) *
+                          oor0[n - nmin];
       }
-      #pragma omp simd
+#pragma omp simd
       for (int n = nmin; n < nmax; ++n) {
         //now r0 depends on ipt and phi as well
-        tmp[n-nmin] = dadipt[n-nmin] * ipt[n-nmin];
+        tmp[n - nmin] = dadipt[n - nmin] * ipt[n - nmin];
       }
-      #pragma omp simd
+#pragma omp simd
       for (int n = nmin; n < nmax; ++n) {
-        dDdipt[n-nmin] -=
-            k[n-nmin] *
-            (x[n-nmin] * (pxca[n-nmin] * tmp[n-nmin] - pysa[n-nmin] * tmp[n-nmin] - pyca[n-nmin] - pxsa[n-nmin] + pyin[n-nmin]) + y[n-nmin] * (pyca[n-nmin] * tmp[n-nmin] + pxsa[n-nmin] * tmp[n-nmin] - pysa[n-nmin] + pxca[n-nmin] - pxin[n-nmin])) *
-            pt[n-nmin] * oor0[n-nmin];
+        dDdipt[n - nmin] -= k[n - nmin] *
+                            (x[n - nmin] * (pxca[n - nmin] * tmp[n - nmin] - pysa[n - nmin] * tmp[n - nmin] -
+                                            pyca[n - nmin] - pxsa[n - nmin] + pyin[n - nmin]) +
+                             y[n - nmin] * (pyca[n - nmin] * tmp[n - nmin] + pxsa[n - nmin] * tmp[n - nmin] -
+                                            pysa[n - nmin] + pxca[n - nmin] - pxin[n - nmin])) *
+                            pt[n - nmin] * oor0[n - nmin];
       }
-      #pragma omp simd
+#pragma omp simd
       for (int n = nmin; n < nmax; ++n) {
-        dDdphi[n-nmin] += k[n-nmin] * (x[n-nmin] * (pysa[n-nmin] - pxin[n-nmin] + pxca[n-nmin]) - y[n-nmin] * (pxsa[n-nmin] - pyin[n-nmin] + pyca[n-nmin])) * oor0[n-nmin];
+        dDdphi[n - nmin] += k[n - nmin] *
+                            (x[n - nmin] * (pysa[n - nmin] - pxin[n - nmin] + pxca[n - nmin]) -
+                             y[n - nmin] * (pxsa[n - nmin] - pyin[n - nmin] + pyca[n - nmin])) *
+                            oor0[n - nmin];
       }
     }
 
-    #pragma omp simd
+#pragma omp simd
     for (int n = nmin; n < nmax; ++n) {
       //update parameters
-      outPar(n, 0, 0) = outPar(n, 0, 0) + 2.f * k[n-nmin] * sinah[n-nmin] * (pxin[n-nmin] * cosah[n-nmin] - pyin[n-nmin] * sinah[n-nmin]);
-      outPar(n, 1, 0) = outPar(n, 1, 0) + 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];
+      outPar(n, 0, 0) = outPar(n, 0, 0) + 2.f * k[n - nmin] * sinah[n - nmin] *
+                                              (pxin[n - nmin] * cosah[n - nmin] - pyin[n - nmin] * sinah[n - nmin]);
+      outPar(n, 1, 0) = outPar(n, 1, 0) + 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];
 
       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
+                "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 alpha[nmax-nmin];
-  float dadphi[nmax-nmin];  
+  float alpha[nmax - nmin];
+  float dadphi[nmax - nmin];
 
-  #pragma omp simd
+#pragma omp simd
   for (int n = nmin; n < nmax; ++n) {
-    alpha[n-nmin] = D[n-nmin] * ipt[n-nmin] * kinv[n-nmin];
-    dadx[n-nmin] = dDdx[n-nmin] * ipt[n-nmin] * kinv[n-nmin];
-    dady[n-nmin] = dDdy[n-nmin] * ipt[n-nmin] * kinv[n-nmin];
-    dadipt[n-nmin] = (ipt[n-nmin] * dDdipt[n-nmin] + D[n-nmin]) * kinv[n-nmin];
-    dadphi[n-nmin] = dDdphi[n-nmin] * ipt[n-nmin] * kinv[n-nmin];
+    alpha[n - nmin] = D[n - nmin] * ipt[n - nmin] * kinv[n - nmin];
+    dadx[n - nmin] = dDdx[n - nmin] * ipt[n - nmin] * kinv[n - nmin];
+    dady[n - nmin] = dDdy[n - nmin] * ipt[n - nmin] * kinv[n - nmin];
+    dadipt[n - nmin] = (ipt[n - nmin] * dDdipt[n - nmin] + D[n - nmin]) * kinv[n - nmin];
+    dadphi[n - nmin] = dDdphi[n - nmin] * ipt[n - nmin] * kinv[n - nmin];
   }
 
   if constexpr (Config::useTrigApprox) {
-    #pragma omp simd
+#pragma omp simd
     for (int n = nmin; n < nmax; ++n) {
-      sincos4(alpha[n-nmin], sina[n-nmin], cosa[n-nmin]);
+      sincos4(alpha[n - nmin], sina[n - nmin], cosa[n - nmin]);
     }
   } else {
-    #pragma omp simd
+#pragma omp simd
     for (int n = nmin; n < nmax; ++n) {
-      cosa[n-nmin] = std::cos(alpha[n-nmin]);
-      sina[n-nmin] = std::sin(alpha[n-nmin]);
+      cosa[n - nmin] = std::cos(alpha[n - nmin]);
+      sina[n - nmin] = std::sin(alpha[n - nmin]);
     }
   }
-  #pragma omp simd
+#pragma omp simd
   for (int n = nmin; n < nmax; ++n) {
-    errorProp(n, 0, 0) = 1.f + k[n-nmin] * dadx[n-nmin] * (cosPorT[n-nmin] * cosa[n-nmin] - sinPorT[n-nmin] * sina[n-nmin]) * pt[n-nmin];
-    errorProp(n, 0, 1) = k[n-nmin] * dady[n-nmin] * (cosPorT[n-nmin] * cosa[n-nmin] - sinPorT[n-nmin] * sina[n-nmin]) * pt[n-nmin];
+    errorProp(n, 0, 0) = 1.f + k[n - nmin] * dadx[n - nmin] *
+                                   (cosPorT[n - nmin] * cosa[n - nmin] - sinPorT[n - nmin] * sina[n - nmin]) *
+                                   pt[n - nmin];
+    errorProp(n, 0, 1) = k[n - nmin] * dady[n - nmin] *
+                         (cosPorT[n - nmin] * cosa[n - nmin] - sinPorT[n - nmin] * sina[n - nmin]) * pt[n - nmin];
     errorProp(n, 0, 2) = 0.f;
     errorProp(n, 0, 3) =
-        k[n-nmin] * (cosPorT[n-nmin] * (ipt[n-nmin] * dadipt[n-nmin] * cosa[n-nmin] - sina[n-nmin]) + sinPorT[n-nmin] * ((1.f - cosa[n-nmin]) - ipt[n-nmin] * dadipt[n-nmin] * sina[n-nmin])) * pt[n-nmin] * pt[n-nmin];
+        k[n - nmin] *
+        (cosPorT[n - nmin] * (ipt[n - nmin] * dadipt[n - nmin] * cosa[n - nmin] - sina[n - nmin]) +
+         sinPorT[n - nmin] * ((1.f - cosa[n - nmin]) - ipt[n - nmin] * dadipt[n - nmin] * sina[n - nmin])) *
+        pt[n - nmin] * pt[n - nmin];
     errorProp(n, 0, 4) =
-        k[n-nmin] * (cosPorT[n-nmin] * dadphi[n-nmin] * cosa[n-nmin] - sinPorT[n-nmin] * dadphi[n-nmin] * sina[n-nmin] - sinPorT[n-nmin] * sina[n-nmin] + cosPorT[n-nmin] * cosa[n-nmin] - cosPorT[n-nmin]) * pt[n-nmin];
+        k[n - nmin] *
+        (cosPorT[n - nmin] * dadphi[n - nmin] * cosa[n - nmin] - sinPorT[n - nmin] * dadphi[n - nmin] * sina[n - nmin] -
+         sinPorT[n - nmin] * sina[n - nmin] + cosPorT[n - nmin] * cosa[n - nmin] - cosPorT[n - nmin]) *
+        pt[n - nmin];
     errorProp(n, 0, 5) = 0.f;
-
   }
 
-  #pragma omp simd
+#pragma omp simd
   for (int n = nmin; n < nmax; ++n) {
-    errorProp(n, 1, 0) = k[n-nmin] * dadx[n-nmin] * (sinPorT[n-nmin] * cosa[n-nmin] + cosPorT[n-nmin] * sina[n-nmin]) * pt[n-nmin];
-    errorProp(n, 1, 1) = 1.f + k[n-nmin] * dady[n-nmin] * (sinPorT[n-nmin] * cosa[n-nmin] + cosPorT[n-nmin] * sina[n-nmin]) * pt[n-nmin];
+    errorProp(n, 1, 0) = k[n - nmin] * dadx[n - nmin] *
+                         (sinPorT[n - nmin] * cosa[n - nmin] + cosPorT[n - nmin] * sina[n - nmin]) * pt[n - nmin];
+    errorProp(n, 1, 1) = 1.f + k[n - nmin] * dady[n - nmin] *
+                                   (sinPorT[n - nmin] * cosa[n - nmin] + cosPorT[n - nmin] * sina[n - nmin]) *
+                                   pt[n - nmin];
     errorProp(n, 1, 2) = 0.f;
     errorProp(n, 1, 3) =
-        k[n-nmin] * (sinPorT[n-nmin] * (ipt[n-nmin] * dadipt[n-nmin] * cosa[n-nmin] - sina[n-nmin]) + cosPorT[n-nmin] * (ipt[n-nmin] * dadipt[n-nmin] * sina[n-nmin] - (1.f - cosa[n-nmin]))) * pt[n-nmin] * pt[n-nmin];
+        k[n - nmin] *
+        (sinPorT[n - nmin] * (ipt[n - nmin] * dadipt[n - nmin] * cosa[n - nmin] - sina[n - nmin]) +
+         cosPorT[n - nmin] * (ipt[n - nmin] * dadipt[n - nmin] * sina[n - nmin] - (1.f - cosa[n - nmin]))) *
+        pt[n - nmin] * pt[n - nmin];
     errorProp(n, 1, 4) =
-        k[n-nmin] * (sinPorT[n-nmin] * dadphi[n-nmin] * cosa[n-nmin] + cosPorT[n-nmin] * dadphi[n-nmin] * sina[n-nmin] + sinPorT[n-nmin] * cosa[n-nmin] + cosPorT[n-nmin] * sina[n-nmin] - sinPorT[n-nmin]) * pt[n-nmin];
+        k[n - nmin] *
+        (sinPorT[n - nmin] * dadphi[n - nmin] * cosa[n - nmin] + cosPorT[n - nmin] * dadphi[n - nmin] * sina[n - nmin] +
+         sinPorT[n - nmin] * cosa[n - nmin] + cosPorT[n - nmin] * sina[n - nmin] - sinPorT[n - nmin]) *
+        pt[n - nmin];
     errorProp(n, 1, 5) = 0.f;
   }
 
-  #pragma omp simd
+#pragma omp simd
   for (int n = nmin; n < nmax; ++n) {
     //no trig approx here, theta can be large
-    cosPorT[n-nmin] = std::cos(theta[n-nmin]);
+    cosPorT[n - nmin] = std::cos(theta[n - nmin]);
   }
 
-  #pragma omp simd
+#pragma omp simd
   for (int n = nmin; n < nmax; ++n) {
-    sinPorT[n-nmin] = std::sin(theta[n-nmin]);
+    sinPorT[n - nmin] = std::sin(theta[n - nmin]);
   }
 
-  #pragma omp simd
+#pragma omp simd
   for (int n = nmin; n < nmax; ++n) {
     //redefine sinPorT as 1./sinPorT to reduce the number of temporaries
-    sinPorT[n-nmin] = 1.f / sinPorT[n-nmin];
+    sinPorT[n - nmin] = 1.f / sinPorT[n - nmin];
   }
 
-  #pragma omp simd
+#pragma omp simd
   for (int n = nmin; n < nmax; ++n) {
-    outPar(n, 2, 0) = inPar(n, 2, 0) + k[n-nmin] * alpha[n-nmin] * cosPorT[n-nmin] * pt[n-nmin] * sinPorT[n-nmin];
-    errorProp(n, 2, 0) = k[n-nmin] * cosPorT[n-nmin] * dadx[n-nmin] * pt[n-nmin] * sinPorT[n-nmin];
-    errorProp(n, 2, 1) = k[n-nmin] * cosPorT[n-nmin] * dady[n-nmin] * pt[n-nmin] * sinPorT[n-nmin];
+    outPar(n, 2, 0) =
+        inPar(n, 2, 0) + k[n - nmin] * alpha[n - nmin] * cosPorT[n - nmin] * pt[n - nmin] * sinPorT[n - nmin];
+    errorProp(n, 2, 0) = k[n - nmin] * cosPorT[n - nmin] * dadx[n - nmin] * pt[n - nmin] * sinPorT[n - nmin];
+    errorProp(n, 2, 1) = k[n - nmin] * cosPorT[n - nmin] * dady[n - nmin] * pt[n - nmin] * sinPorT[n - nmin];
     errorProp(n, 2, 2) = 1.f;
-    errorProp(n, 2, 3) = k[n-nmin] * cosPorT[n-nmin] * (ipt[n-nmin] * dadipt[n-nmin] - alpha[n-nmin]) * pt[n-nmin] * pt[n-nmin] * sinPorT[n-nmin];
-    errorProp(n, 2, 4) = k[n-nmin] * dadphi[n-nmin] * cosPorT[n-nmin] * pt[n-nmin] * sinPorT[n-nmin];
-    errorProp(n, 2, 5) = -k[n-nmin] * alpha[n-nmin] * pt[n-nmin] * sinPorT[n-nmin] * sinPorT[n-nmin];
+    errorProp(n, 2, 3) = k[n - nmin] * cosPorT[n - nmin] * (ipt[n - nmin] * dadipt[n - nmin] - alpha[n - nmin]) *
+                         pt[n - nmin] * pt[n - nmin] * sinPorT[n - nmin];
+    errorProp(n, 2, 4) = k[n - nmin] * dadphi[n - nmin] * cosPorT[n - nmin] * pt[n - nmin] * sinPorT[n - nmin];
+    errorProp(n, 2, 5) = -k[n - nmin] * alpha[n - nmin] * pt[n - nmin] * sinPorT[n - nmin] * sinPorT[n - nmin];
   }
 
-  #pragma omp simd
+#pragma omp simd
   for (int n = nmin; n < nmax; ++n) {
-    outPar(n, 3, 0) = ipt[n-nmin];
+    outPar(n, 3, 0) = ipt[n - nmin];
     errorProp(n, 3, 0) = 0.f;
     errorProp(n, 3, 1) = 0.f;
     errorProp(n, 3, 2) = 0.f;
@@ -381,20 +413,20 @@ static inline void helixAtRFromIterativeCCS_impl(const Tf& __restrict__ inPar,
     errorProp(n, 3, 5) = 0.f;
   }
 
-  #pragma omp simd
+#pragma omp simd
   for (int n = nmin; n < nmax; ++n) {
-    outPar(n, 4, 0) = inPar(n, 4, 0) + alpha[n-nmin];
-    errorProp(n, 4, 0) = dadx[n-nmin];
-    errorProp(n, 4, 1) = dady[n-nmin];
+    outPar(n, 4, 0) = inPar(n, 4, 0) + alpha[n - nmin];
+    errorProp(n, 4, 0) = dadx[n - nmin];
+    errorProp(n, 4, 1) = dady[n - nmin];
     errorProp(n, 4, 2) = 0.f;
-    errorProp(n, 4, 3) = dadipt[n-nmin];
-    errorProp(n, 4, 4) = 1.f + dadphi[n-nmin];
+    errorProp(n, 4, 3) = dadipt[n - nmin];
+    errorProp(n, 4, 4) = 1.f + dadphi[n - nmin];
     errorProp(n, 4, 5) = 0.f;
   }
 
-  #pragma omp simd
+#pragma omp simd
   for (int n = nmin; n < nmax; ++n) {
-    outPar(n, 5, 0) = theta[n-nmin];
+    outPar(n, 5, 0) = theta[n - nmin];
     errorProp(n, 5, 0) = 0.f;
     errorProp(n, 5, 1) = 0.f;
     errorProp(n, 5, 2) = 0.f;
@@ -403,7 +435,7 @@ static inline void helixAtRFromIterativeCCS_impl(const Tf& __restrict__ inPar,
     errorProp(n, 5, 5) = 1.f;
   }
 
-  #pragma omp simd
+#pragma omp simd
   for (int n = nmin; n < nmax; ++n) {
     dprint_np(n,
               "propagation end, dump parameters"