MkFitCore/src/PropagationMPlex.icc

0001 ///////////////////////////////////////////////////////////////////////////////
0002 /// helixAtRFromIterativeCCS_impl
0003 ///////////////////////////////////////////////////////////////////////////////
0004
0005 template <typename Tf, typename Ti, typename TfLL1, typename Tf11, typename TfLLL>
0006 static inline void helixAtRFromIterativeCCS_impl(const Tf& __restrict__ inPar,
0007                                                  const Ti& __restrict__ inChg,
0008                                                  const Tf11& __restrict__ msRad,
0009                                                  TfLL1& __restrict__ outPar,
0010                                                  TfLLL& __restrict__ errorProp,
0011                                                  Ti& __restrict__ outFailFlag,  // expected to be initialized to 0
0012                                                  const int nmin,
0013                                                  const int nmax,
0014                                                  const int N_proc,
0015                                                  const PropagationFlags pf) {
0016   // bool debug = true;
0017
0018 #pragma omp simd
0019   for (int n = nmin; n < nmax; ++n) {
0020     //initialize erroProp to identity matrix
0021     errorProp(n, 0, 0) = 1.f;
0022     errorProp(n, 1, 1) = 1.f;
0023     errorProp(n, 2, 2) = 1.f;
0024     errorProp(n, 3, 3) = 1.f;
0025     errorProp(n, 4, 4) = 1.f;
0026     errorProp(n, 5, 5) = 1.f;
0027   }
0028   float r0[nmax - nmin];
0029 #pragma omp simd
0030   for (int n = nmin; n < nmax; ++n) {
0031     //initialize erroProp to identity matrix
0032     r0[n - nmin] = hipo(inPar(n, 0, 0), inPar(n, 1, 0));
0033   }
0034   float k[nmax - nmin];
0035   if (pf.use_param_b_field) {
0036 #pragma omp simd
0037     for (int n = nmin; n < nmax; ++n) {
0038       k[n - nmin] = inChg(n, 0, 0) * 100.f / (-Const::sol * Config::bFieldFromZR(inPar(n, 2, 0), r0[n - nmin]));
0039     }
0040   } else {
0041 #pragma omp simd
0042     for (int n = nmin; n < nmax; ++n) {
0043       k[n - nmin] = inChg(n, 0, 0) * 100.f / (-Const::sol * Config::Bfield);
0044     }
0045   }
0046   float r[nmax - nmin];
0047 #pragma omp simd
0048   for (int n = nmin; n < nmax; ++n) {
0049     r[n - nmin] = msRad(n, 0, 0);
0050   }
0051   float xin[nmax - nmin];
0052   float yin[nmax - nmin];
0053   float ipt[nmax - nmin];
0054   float phiin[nmax - nmin];
0055   float theta[nmax - nmin];
0056 #pragma omp simd
0057   for (int n = nmin; n < nmax; ++n) {
0058     // if (std::abs(r-r0)<0.0001f) {
0059     //  dprint("distance less than 1mum, skip");
0060     //  continue;
0061     // }
0062
0063     xin[n - nmin] = inPar(n, 0, 0);
0064     yin[n - nmin] = inPar(n, 1, 0);
0065     ipt[n - nmin] = inPar(n, 3, 0);
0066     phiin[n - nmin] = inPar(n, 4, 0);
0067     theta[n - nmin] = inPar(n, 5, 0);
0068
0069     //dprint(std::endl);
0070   }
0071
0072   for (int n = nmin; n < nmax; ++n) {
0073     dprint_np(n,
0074               "input parameters"
0075                   << " inPar(n, 0, 0)=" << std::setprecision(9) << inPar(n, 0, 0) << " inPar(n, 1, 0)="
0076                   << std::setprecision(9) << inPar(n, 1, 0) << " inPar(n, 2, 0)=" << std::setprecision(9)
0077                   << inPar(n, 2, 0) << " inPar(n, 3, 0)=" << std::setprecision(9) << inPar(n, 3, 0)
0078                   << " inPar(n, 4, 0)=" << std::setprecision(9) << inPar(n, 4, 0)
0079                   << " inPar(n, 5, 0)=" << std::setprecision(9) << inPar(n, 5, 0));
0080   }
0081
0082   float kinv[nmax - nmin];
0083   float pt[nmax - nmin];
0084 #pragma omp simd
0085   for (int n = nmin; n < nmax; ++n) {
0086     kinv[n - nmin] = 1.f / k[n - nmin];
0087     pt[n - nmin] = 1.f / ipt[n - nmin];
0088   }
0089   float D[nmax - nmin];
0090   float cosa[nmax - nmin];
0091   float sina[nmax - nmin];
0092   float cosah[nmax - nmin];
0093   float sinah[nmax - nmin];
0094   float id[nmax - nmin];
0095
0096 #pragma omp simd
0097   for (int n = nmin; n < nmax; ++n) {
0098     D[n - nmin] = 0.;
0099   }
0100
0101   //no trig approx here, phi can be large
0102   float cosPorT[nmax - nmin];
0103   float sinPorT[nmax - nmin];
0104 #pragma omp simd
0105   for (int n = nmin; n < nmax; ++n) {
0106     cosPorT[n - nmin] = std::cos(phiin[n - nmin]);
0107   }
0108 #pragma omp simd
0109   for (int n = nmin; n < nmax; ++n) {
0110     sinPorT[n - nmin] = std::sin(phiin[n - nmin]);
0111   }
0112
0113   float pxin[nmax - nmin];
0114   float pyin[nmax - nmin];
0115 #pragma omp simd
0116   for (int n = nmin; n < nmax; ++n) {
0117     pxin[n - nmin] = cosPorT[n - nmin] * pt[n - nmin];
0118     pyin[n - nmin] = sinPorT[n - nmin] * pt[n - nmin];
0119   }
0120
0121   for (int n = nmin; n < nmax; ++n) {
0122     dprint_np(n,
0123               "k=" << std::setprecision(9) << k[n - nmin] << " pxin=" << std::setprecision(9) << pxin[n - nmin]
0124                    << " pyin=" << std::setprecision(9) << pyin[n - nmin] << " cosPorT=" << std::setprecision(9)
0125                    << cosPorT[n - nmin] << " sinPorT=" << std::setprecision(9) << sinPorT[n - nmin]
0126                    << " pt=" << std::setprecision(9) << pt[n - nmin]);
0127   }
0128
0129   float dDdx[nmax - nmin];
0130   float dDdy[nmax - nmin];
0131   float dDdipt[nmax - nmin];
0132   float dDdphi[nmax - nmin];
0133
0134 #pragma omp simd
0135   for (int n = nmin; n < nmax; ++n) {
0136     dDdipt[n - nmin] = 0.;
0137     dDdphi[n - nmin] = 0.;
0138   }
0139 #pragma omp simd
0140   for (int n = nmin; n < nmax; ++n) {
0141     //derivatives initialized to value for first iteration, i.e. distance = r-r0in
0142     dDdx[n - nmin] = r0[n - nmin] > 0.f ? -xin[n - nmin] / r0[n - nmin] : 0.f;
0143   }
0144
0145 #pragma omp simd
0146   for (int n = nmin; n < nmax; ++n) {
0147     dDdy[n - nmin] = r0[n - nmin] > 0.f ? -yin[n - nmin] / r0[n - nmin] : 0.f;
0148   }
0149
0150   float oodotp[nmax - nmin];
0151   float x[nmax - nmin];
0152   float y[nmax - nmin];
0153   float oor0[nmax - nmin];
0154   float dadipt[nmax - nmin];
0155   float dadx[nmax - nmin];
0156   float dady[nmax - nmin];
0157   float pxca[nmax - nmin];
0158   float pxsa[nmax - nmin];
0159   float pyca[nmax - nmin];
0160   float pysa[nmax - nmin];
0161   float tmp[nmax - nmin];
0162   float tmpx[nmax - nmin];
0163   float tmpy[nmax - nmin];
0164   float pxinold[nmax - nmin];
0165
0166   CMS_UNROLL_LOOP_COUNT(Config::Niter)
0167   for (int i = 0; i < Config::Niter; ++i) {
0168 #pragma omp simd
0169     for (int n = nmin; n < nmax; ++n) {
0170       //compute distance and path for the current iteration
0171       r0[n - nmin] = hipo(outPar(n, 0, 0), outPar(n, 1, 0));
0172     }
0173
0174     // Use one over dot produce of transverse momentum and radial
0175     // direction to scale the step. Propagation is prevented from reaching
0176     // too close to the apex (dotp > 0.2).
0177     // - Can / should we come up with a better approximation?
0178     // - Can / should take +/- curvature into account?
0179
0180 #pragma omp simd
0181     for (int n = nmin; n < nmax; ++n) {
0182       oodotp[n - nmin] =
0183           r0[n - nmin] * pt[n - nmin] / (pxin[n - nmin] * outPar(n, 0, 0) + pyin[n - nmin] * outPar(n, 1, 0));
0184     }
0185
0186 #pragma omp simd
0187     for (int n = nmin; n < nmax; ++n) {
0188       if (oodotp[n - nmin] > 5.0f || oodotp[n - nmin] < 0)  // 0.2 is 78.5 deg
0189       {
0190         outFailFlag(n, 0, 0) = 1;
0191       }
0192     }
0193
0194 #pragma omp simd
0195     for (int n = nmin; n < nmax; ++n) {
0196       // Can we come up with a better approximation?
0197       // Should take +/- curvature into account.
0198       id[n - nmin] =
0199           (oodotp[n - nmin] > 5.0f || oodotp[n - nmin] < 0) ? 0.0f : (r[n - nmin] - r0[n - nmin]) * oodotp[n - nmin];
0200     }
0201
0202 #pragma omp simd
0203     for (int n = nmin; n < nmax; ++n) {
0204       D[n - nmin] += id[n - nmin];
0205     }
0206
0207     if constexpr (Config::useTrigApprox) {
0208 #if !defined(__INTEL_COMPILER)
0209 #pragma omp simd
0210 #endif
0211       for (int n = nmin; n < nmax; ++n) {
0212         sincos4(id[n - nmin] * ipt[n - nmin] * kinv[n - nmin] * 0.5f, sinah[n - nmin], cosah[n - nmin]);
0213       }
0214     } else {
0215 #if !defined(__INTEL_COMPILER)
0216 #pragma omp simd
0217 #endif
0218       for (int n = nmin; n < nmax; ++n) {
0219         cosah[n - nmin] = std::cos(id[n - nmin] * ipt[n - nmin] * kinv[n - nmin] * 0.5f);
0220         sinah[n - nmin] = std::sin(id[n - nmin] * ipt[n - nmin] * kinv[n - nmin] * 0.5f);
0221       }
0222     }
0223
0224 #pragma omp simd
0225     for (int n = nmin; n < nmax; ++n) {
0226       cosa[n - nmin] = 1.f - 2.f * sinah[n - nmin] * sinah[n - nmin];
0227       sina[n - nmin] = 2.f * sinah[n - nmin] * cosah[n - nmin];
0228     }
0229
0230     for (int n = nmin; n < nmax; ++n) {
0231       dprint_np(n,
0232                 "Attempt propagation from r="
0233                     << r0[n - nmin] << " to r=" << r[n - nmin] << std::endl
0234                     << "   x=" << xin[n - nmin] << " y=" << yin[n - nmin] << " z=" << inPar(n, 2, 0)
0235                     << " px=" << pxin[n - nmin] << " py=" << pyin[n - nmin]
0236                     << " pz=" << pt[n - nmin] * std::tan(theta[n - nmin]) << " q=" << inChg(n, 0, 0) << std::endl
0237                     << "   r=" << std::setprecision(9) << r[n - nmin] << " r0=" << std::setprecision(9) << r0[n - nmin]
0238                     << " id=" << std::setprecision(9) << id[n - nmin] << " dr=" << std::setprecision(9)
0239                     << r[n - nmin] - r0[n - nmin] << " cosa=" << cosa[n - nmin] << " sina=" << sina[n - nmin]);
0240     }
0241
0242     //update derivatives on total distance
0243     if (i + 1 != Config::Niter) {
0244 #pragma omp simd
0245       for (int n = nmin; n < nmax; ++n) {
0246         x[n - nmin] = outPar(n, 0, 0);
0247         y[n - nmin] = outPar(n, 1, 0);
0248       }
0249 #pragma omp simd
0250       for (int n = nmin; n < nmax; ++n) {
0251         oor0[n - nmin] =
0252             (r0[n - nmin] > 0.f && std::abs(r[n - nmin] - r0[n - nmin]) < 0.0001f) ? 1.f / r0[n - nmin] : 0.f;
0253       }
0254 #pragma omp simd
0255       for (int n = nmin; n < nmax; ++n) {
0256         dadipt[n - nmin] = id[n - nmin] * kinv[n - nmin];
0257         dadx[n - nmin] = -x[n - nmin] * ipt[n - nmin] * kinv[n - nmin] * oor0[n - nmin];
0258         dady[n - nmin] = -y[n - nmin] * ipt[n - nmin] * kinv[n - nmin] * oor0[n - nmin];
0259         pxca[n - nmin] = pxin[n - nmin] * cosa[n - nmin];
0260         pxsa[n - nmin] = pxin[n - nmin] * sina[n - nmin];
0261         pyca[n - nmin] = pyin[n - nmin] * cosa[n - nmin];
0262         pysa[n - nmin] = pyin[n - nmin] * sina[n - nmin];
0263         tmpx[n - nmin] = k[n - nmin] * dadx[n - nmin];
0264       }
0265
0266 #pragma omp simd
0267       for (int n = nmin; n < nmax; ++n) {
0268         dDdx[n - nmin] -= (x[n - nmin] * (1.f + tmpx[n - nmin] * (pxca[n - nmin] - pysa[n - nmin])) +
0269                            y[n - nmin] * tmpx[n - nmin] * (pyca[n - nmin] + pxsa[n - nmin])) *
0270                           oor0[n - nmin];
0271       }
0272
0273 #pragma omp simd
0274       for (int n = nmin; n < nmax; ++n) {
0275         tmpy[n - nmin] = k[n - nmin] * dady[n - nmin];
0276       }
0277 #pragma omp simd
0278       for (int n = nmin; n < nmax; ++n) {
0279         dDdy[n - nmin] -= (x[n - nmin] * tmpy[n - nmin] * (pxca[n - nmin] - pysa[n - nmin]) +
0280                            y[n - nmin] * (1.f + tmpy[n - nmin] * (pyca[n - nmin] + pxsa[n - nmin]))) *
0281                           oor0[n - nmin];
0282       }
0283 #pragma omp simd
0284       for (int n = nmin; n < nmax; ++n) {
0285         //now r0 depends on ipt and phi as well
0286         tmp[n - nmin] = dadipt[n - nmin] * ipt[n - nmin];
0287       }
0288 #pragma omp simd
0289       for (int n = nmin; n < nmax; ++n) {
0290         dDdipt[n - nmin] -= k[n - nmin] *
0291                             (x[n - nmin] * (pxca[n - nmin] * tmp[n - nmin] - pysa[n - nmin] * tmp[n - nmin] -
0292                                             pyca[n - nmin] - pxsa[n - nmin] + pyin[n - nmin]) +
0293                              y[n - nmin] * (pyca[n - nmin] * tmp[n - nmin] + pxsa[n - nmin] * tmp[n - nmin] -
0294                                             pysa[n - nmin] + pxca[n - nmin] - pxin[n - nmin])) *
0295                             pt[n - nmin] * oor0[n - nmin];
0296       }
0297 #pragma omp simd
0298       for (int n = nmin; n < nmax; ++n) {
0299         dDdphi[n - nmin] += k[n - nmin] *
0300                             (x[n - nmin] * (pysa[n - nmin] - pxin[n - nmin] + pxca[n - nmin]) -
0301                              y[n - nmin] * (pxsa[n - nmin] - pyin[n - nmin] + pyca[n - nmin])) *
0302                             oor0[n - nmin];
0303       }
0304     }
0305
0306 #pragma omp simd
0307     for (int n = nmin; n < nmax; ++n) {
0308       //update parameters
0309       outPar(n, 0, 0) = outPar(n, 0, 0) + 2.f * k[n - nmin] * sinah[n - nmin] *
0310                                               (pxin[n - nmin] * cosah[n - nmin] - pyin[n - nmin] * sinah[n - nmin]);
0311       outPar(n, 1, 0) = outPar(n, 1, 0) + 2.f * k[n - nmin] * sinah[n - nmin] *
0312                                               (pyin[n - nmin] * cosah[n - nmin] + pxin[n - nmin] * sinah[n - nmin]);
0313       pxinold[n - nmin] = pxin[n - nmin];  //copy before overwriting
0314       pxin[n - nmin] = pxin[n - nmin] * cosa[n - nmin] - pyin[n - nmin] * sina[n - nmin];
0315       pyin[n - nmin] = pyin[n - nmin] * cosa[n - nmin] + pxinold[n - nmin] * sina[n - nmin];
0316     }
0317     for (int n = nmin; n < nmax; ++n) {
0318       dprint_np(n,
0319                 "outPar(n, 0, 0)=" << outPar(n, 0, 0) << " outPar(n, 1, 0)=" << outPar(n, 1, 0)
0320                                    << " pxin=" << pxin[n - nmin] << " pyin=" << pyin[n - nmin]);
0321     }
0322   }  // iteration loop
0323
0324   float alpha[nmax - nmin];
0325   float dadphi[nmax - nmin];
0326
0327 #pragma omp simd
0328   for (int n = nmin; n < nmax; ++n) {
0329     alpha[n - nmin] = D[n - nmin] * ipt[n - nmin] * kinv[n - nmin];
0330     dadx[n - nmin] = dDdx[n - nmin] * ipt[n - nmin] * kinv[n - nmin];
0331     dady[n - nmin] = dDdy[n - nmin] * ipt[n - nmin] * kinv[n - nmin];
0332     dadipt[n - nmin] = (ipt[n - nmin] * dDdipt[n - nmin] + D[n - nmin]) * kinv[n - nmin];
0333     dadphi[n - nmin] = dDdphi[n - nmin] * ipt[n - nmin] * kinv[n - nmin];
0334   }
0335
0336   if constexpr (Config::useTrigApprox) {
0337 #pragma omp simd
0338     for (int n = nmin; n < nmax; ++n) {
0339       sincos4(alpha[n - nmin], sina[n - nmin], cosa[n - nmin]);
0340     }
0341   } else {
0342 #pragma omp simd
0343     for (int n = nmin; n < nmax; ++n) {
0344       cosa[n - nmin] = std::cos(alpha[n - nmin]);
0345       sina[n - nmin] = std::sin(alpha[n - nmin]);
0346     }
0347   }
0348 #pragma omp simd
0349   for (int n = nmin; n < nmax; ++n) {
0350     errorProp(n, 0, 0) = 1.f + k[n - nmin] * dadx[n - nmin] *
0351                                    (cosPorT[n - nmin] * cosa[n - nmin] - sinPorT[n - nmin] * sina[n - nmin]) *
0352                                    pt[n - nmin];
0353     errorProp(n, 0, 1) = k[n - nmin] * dady[n - nmin] *
0354                          (cosPorT[n - nmin] * cosa[n - nmin] - sinPorT[n - nmin] * sina[n - nmin]) * pt[n - nmin];
0355     errorProp(n, 0, 2) = 0.f;
0356     errorProp(n, 0, 3) =
0357         k[n - nmin] *
0358         (cosPorT[n - nmin] * (ipt[n - nmin] * dadipt[n - nmin] * cosa[n - nmin] - sina[n - nmin]) +
0359          sinPorT[n - nmin] * ((1.f - cosa[n - nmin]) - ipt[n - nmin] * dadipt[n - nmin] * sina[n - nmin])) *
0360         pt[n - nmin] * pt[n - nmin];
0361     errorProp(n, 0, 4) =
0362         k[n - nmin] *
0363         (cosPorT[n - nmin] * dadphi[n - nmin] * cosa[n - nmin] - sinPorT[n - nmin] * dadphi[n - nmin] * sina[n - nmin] -
0364          sinPorT[n - nmin] * sina[n - nmin] + cosPorT[n - nmin] * cosa[n - nmin] - cosPorT[n - nmin]) *
0365         pt[n - nmin];
0366     errorProp(n, 0, 5) = 0.f;
0367   }
0368
0369 #pragma omp simd
0370   for (int n = nmin; n < nmax; ++n) {
0371     errorProp(n, 1, 0) = k[n - nmin] * dadx[n - nmin] *
0372                          (sinPorT[n - nmin] * cosa[n - nmin] + cosPorT[n - nmin] * sina[n - nmin]) * pt[n - nmin];
0373     errorProp(n, 1, 1) = 1.f + k[n - nmin] * dady[n - nmin] *
0374                                    (sinPorT[n - nmin] * cosa[n - nmin] + cosPorT[n - nmin] * sina[n - nmin]) *
0375                                    pt[n - nmin];
0376     errorProp(n, 1, 2) = 0.f;
0377     errorProp(n, 1, 3) =
0378         k[n - nmin] *
0379         (sinPorT[n - nmin] * (ipt[n - nmin] * dadipt[n - nmin] * cosa[n - nmin] - sina[n - nmin]) +
0380          cosPorT[n - nmin] * (ipt[n - nmin] * dadipt[n - nmin] * sina[n - nmin] - (1.f - cosa[n - nmin]))) *
0381         pt[n - nmin] * pt[n - nmin];
0382     errorProp(n, 1, 4) =
0383         k[n - nmin] *
0384         (sinPorT[n - nmin] * dadphi[n - nmin] * cosa[n - nmin] + cosPorT[n - nmin] * dadphi[n - nmin] * sina[n - nmin] +
0385          sinPorT[n - nmin] * cosa[n - nmin] + cosPorT[n - nmin] * sina[n - nmin] - sinPorT[n - nmin]) *
0386         pt[n - nmin];
0387     errorProp(n, 1, 5) = 0.f;
0388   }
0389
0390 #pragma omp simd
0391   for (int n = nmin; n < nmax; ++n) {
0392     //no trig approx here, theta can be large
0393     cosPorT[n - nmin] = std::cos(theta[n - nmin]);
0394   }
0395
0396 #pragma omp simd
0397   for (int n = nmin; n < nmax; ++n) {
0398     sinPorT[n - nmin] = std::sin(theta[n - nmin]);
0399   }
0400
0401 #pragma omp simd
0402   for (int n = nmin; n < nmax; ++n) {
0403     //redefine sinPorT as 1./sinPorT to reduce the number of temporaries
0404     sinPorT[n - nmin] = 1.f / sinPorT[n - nmin];
0405   }
0406
0407 #pragma omp simd
0408   for (int n = nmin; n < nmax; ++n) {
0409     outPar(n, 2, 0) =
0410         inPar(n, 2, 0) + k[n - nmin] * alpha[n - nmin] * cosPorT[n - nmin] * pt[n - nmin] * sinPorT[n - nmin];
0411     errorProp(n, 2, 0) = k[n - nmin] * cosPorT[n - nmin] * dadx[n - nmin] * pt[n - nmin] * sinPorT[n - nmin];
0412     errorProp(n, 2, 1) = k[n - nmin] * cosPorT[n - nmin] * dady[n - nmin] * pt[n - nmin] * sinPorT[n - nmin];
0413     errorProp(n, 2, 2) = 1.f;
0414     errorProp(n, 2, 3) = k[n - nmin] * cosPorT[n - nmin] * (ipt[n - nmin] * dadipt[n - nmin] - alpha[n - nmin]) *
0415                          pt[n - nmin] * pt[n - nmin] * sinPorT[n - nmin];
0416     errorProp(n, 2, 4) = k[n - nmin] * dadphi[n - nmin] * cosPorT[n - nmin] * pt[n - nmin] * sinPorT[n - nmin];
0417     errorProp(n, 2, 5) = -k[n - nmin] * alpha[n - nmin] * pt[n - nmin] * sinPorT[n - nmin] * sinPorT[n - nmin];
0418   }
0419
0420 #pragma omp simd
0421   for (int n = nmin; n < nmax; ++n) {
0422     outPar(n, 3, 0) = ipt[n - nmin];
0423     errorProp(n, 3, 0) = 0.f;
0424     errorProp(n, 3, 1) = 0.f;
0425     errorProp(n, 3, 2) = 0.f;
0426     errorProp(n, 3, 3) = 1.f;
0427     errorProp(n, 3, 4) = 0.f;
0428     errorProp(n, 3, 5) = 0.f;
0429   }
0430
0431 #pragma omp simd
0432   for (int n = nmin; n < nmax; ++n) {
0433     outPar(n, 4, 0) = inPar(n, 4, 0) + alpha[n - nmin];
0434     errorProp(n, 4, 0) = dadx[n - nmin];
0435     errorProp(n, 4, 1) = dady[n - nmin];
0436     errorProp(n, 4, 2) = 0.f;
0437     errorProp(n, 4, 3) = dadipt[n - nmin];
0438     errorProp(n, 4, 4) = 1.f + dadphi[n - nmin];
0439     errorProp(n, 4, 5) = 0.f;
0440   }
0441
0442 #pragma omp simd
0443   for (int n = nmin; n < nmax; ++n) {
0444     outPar(n, 5, 0) = theta[n - nmin];
0445     errorProp(n, 5, 0) = 0.f;
0446     errorProp(n, 5, 1) = 0.f;
0447     errorProp(n, 5, 2) = 0.f;
0448     errorProp(n, 5, 3) = 0.f;
0449     errorProp(n, 5, 4) = 0.f;
0450     errorProp(n, 5, 5) = 1.f;
0451   }
0452
0453   for (int n = nmin; n < nmax; ++n) {
0454     dprint_np(n,
0455               "propagation end, dump parameters"
0456                   << std::endl
0457                   << "   pos = " << outPar(n, 0, 0) << " " << outPar(n, 1, 0) << " " << outPar(n, 2, 0) << "\t\t r="
0458                   << std::sqrt(outPar(n, 0, 0) * outPar(n, 0, 0) + outPar(n, 1, 0) * outPar(n, 1, 0)) << std::endl
0459                   << "   mom = " << std::cos(outPar(n, 4, 0)) / outPar(n, 3, 0) << " "
0460                   << std::sin(outPar(n, 4, 0)) / outPar(n, 3, 0) << " " << 1. / (outPar(n, 3, 0) * tan(outPar(n, 5, 0)))
0461                   << "\t\tpT=" << 1. / std::abs(outPar(n, 3, 0)) << std::endl);
0462   }
0463
0464 #ifdef DEBUG
0465   for (int n = nmin; n < nmax; ++n) {
0466     if (n < N_proc) {
0467       dmutex_guard;
0468       std::cout << n << ": jacobian" << std::endl;
0469       printf("%5f %5f %5f %5f %5f %5f\n",
0470              errorProp(n, 0, 0),
0471              errorProp(n, 0, 1),
0472              errorProp(n, 0, 2),
0473              errorProp(n, 0, 3),
0474              errorProp(n, 0, 4),
0475              errorProp(n, 0, 5));
0476       printf("%5f %5f %5f %5f %5f %5f\n",
0477              errorProp(n, 1, 0),
0478              errorProp(n, 1, 1),
0479              errorProp(n, 1, 2),
0480              errorProp(n, 1, 3),
0481              errorProp(n, 1, 4),
0482              errorProp(n, 1, 5));
0483       printf("%5f %5f %5f %5f %5f %5f\n",
0484              errorProp(n, 2, 0),
0485              errorProp(n, 2, 1),
0486              errorProp(n, 2, 2),
0487              errorProp(n, 2, 3),
0488              errorProp(n, 2, 4),
0489              errorProp(n, 2, 5));
0490       printf("%5f %5f %5f %5f %5f %5f\n",
0491              errorProp(n, 3, 0),
0492              errorProp(n, 3, 1),
0493              errorProp(n, 3, 2),
0494              errorProp(n, 3, 3),
0495              errorProp(n, 3, 4),
0496              errorProp(n, 3, 5));
0497       printf("%5f %5f %5f %5f %5f %5f\n",
0498              errorProp(n, 4, 0),
0499              errorProp(n, 4, 1),
0500              errorProp(n, 4, 2),
0501              errorProp(n, 4, 3),
0502              errorProp(n, 4, 4),
0503              errorProp(n, 4, 5));
0504       printf("%5f %5f %5f %5f %5f %5f\n",
0505              errorProp(n, 5, 0),
0506              errorProp(n, 5, 1),
0507              errorProp(n, 5, 2),
0508              errorProp(n, 5, 3),
0509              errorProp(n, 5, 4),
0510              errorProp(n, 5, 5));
0511       printf("\n");
0512     }
0513   }
0514 #endif
0515 }