Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-05-29 02:55:21 

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 product 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         oodotp[n - nmin] = 0.0f;
0192       } else if (r[n - nmin] - r0[n - nmin] < 0.0f && pt[n - nmin] < 1.0f) {
0193         // Scale down the correction for low-pT ingoing tracks.
0194         oodotp[n - nmin] = 1.0f + (oodotp[n - nmin] - 1.0f) * pt[n - nmin];
0195       }
0196     }
0197 
0198 #pragma omp simd
0199     for (int n = nmin; n < nmax; ++n) {
0200       // Can we come up with a better approximation?
0201       // Should take +/- curvature into account.
0202       id[n - nmin] = (r[n - nmin] - r0[n - nmin]) * oodotp[n - nmin];
0203     }
0204 
0205 #pragma omp simd
0206     for (int n = nmin; n < nmax; ++n) {
0207       D[n - nmin] += id[n - nmin];
0208     }
0209 
0210     if constexpr (Config::useTrigApprox) {
0211 #if !defined(__INTEL_COMPILER)
0212 #pragma omp simd
0213 #endif
0214       for (int n = nmin; n < nmax; ++n) {
0215         sincos4(id[n - nmin] * ipt[n - nmin] * kinv[n - nmin] * 0.5f, sinah[n - nmin], cosah[n - nmin]);
0216       }
0217     } else {
0218 #if !defined(__INTEL_COMPILER)
0219 #pragma omp simd
0220 #endif
0221       for (int n = nmin; n < nmax; ++n) {
0222         cosah[n - nmin] = std::cos(id[n - nmin] * ipt[n - nmin] * kinv[n - nmin] * 0.5f);
0223         sinah[n - nmin] = std::sin(id[n - nmin] * ipt[n - nmin] * kinv[n - nmin] * 0.5f);
0224       }
0225     }
0226 
0227 #pragma omp simd
0228     for (int n = nmin; n < nmax; ++n) {
0229       cosa[n - nmin] = 1.f - 2.f * sinah[n - nmin] * sinah[n - nmin];
0230       sina[n - nmin] = 2.f * sinah[n - nmin] * cosah[n - nmin];
0231     }
0232 
0233     for (int n = nmin; n < nmax; ++n) {
0234       dprint_np(n,
0235                 "Attempt propagation from r="
0236                     << r0[n - nmin] << " to r=" << r[n - nmin] << std::endl
0237                     << "   x=" << xin[n - nmin] << " y=" << yin[n - nmin] << " z=" << inPar(n, 2, 0)
0238                     << " px=" << pxin[n - nmin] << " py=" << pyin[n - nmin]
0239                     << " pz=" << pt[n - nmin] * std::tan(theta[n - nmin]) << " q=" << inChg(n, 0, 0) << std::endl
0240                     << "   r=" << std::setprecision(9) << r[n - nmin] << " r0=" << std::setprecision(9) << r0[n - nmin]
0241                     << " id=" << std::setprecision(9) << id[n - nmin] << " dr=" << std::setprecision(9)
0242                     << r[n - nmin] - r0[n - nmin] << " cosa=" << cosa[n - nmin] << " sina=" << sina[n - nmin]
0243                     << " dir_cos(rad,pT)=" << 1.0f / oodotp[n]);
0244     }
0245 
0246     //update derivatives on total distance
0247     if (i + 1 != Config::Niter) {
0248 #pragma omp simd
0249       for (int n = nmin; n < nmax; ++n) {
0250         x[n - nmin] = outPar(n, 0, 0);
0251         y[n - nmin] = outPar(n, 1, 0);
0252       }
0253 #pragma omp simd
0254       for (int n = nmin; n < nmax; ++n) {
0255         oor0[n - nmin] =
0256             (r0[n - nmin] > 0.f && std::abs(r[n - nmin] - r0[n - nmin]) > 0.0001f) ? 1.f / r0[n - nmin] : 0.f;
0257       }
0258 #pragma omp simd
0259       for (int n = nmin; n < nmax; ++n) {
0260         dadipt[n - nmin] = id[n - nmin] * kinv[n - nmin];
0261         dadx[n - nmin] = -x[n - nmin] * ipt[n - nmin] * kinv[n - nmin] * oor0[n - nmin];
0262         dady[n - nmin] = -y[n - nmin] * ipt[n - nmin] * kinv[n - nmin] * oor0[n - nmin];
0263         pxca[n - nmin] = pxin[n - nmin] * cosa[n - nmin];
0264         pxsa[n - nmin] = pxin[n - nmin] * sina[n - nmin];
0265         pyca[n - nmin] = pyin[n - nmin] * cosa[n - nmin];
0266         pysa[n - nmin] = pyin[n - nmin] * sina[n - nmin];
0267         tmpx[n - nmin] = k[n - nmin] * dadx[n - nmin];
0268       }
0269 
0270 #pragma omp simd
0271       for (int n = nmin; n < nmax; ++n) {
0272         dDdx[n - nmin] -= (x[n - nmin] * (1.f + tmpx[n - nmin] * (pxca[n - nmin] - pysa[n - nmin])) +
0273                            y[n - nmin] * tmpx[n - nmin] * (pyca[n - nmin] + pxsa[n - nmin])) *
0274                           oor0[n - nmin];
0275       }
0276 
0277 #pragma omp simd
0278       for (int n = nmin; n < nmax; ++n) {
0279         tmpy[n - nmin] = k[n - nmin] * dady[n - nmin];
0280       }
0281 #pragma omp simd
0282       for (int n = nmin; n < nmax; ++n) {
0283         dDdy[n - nmin] -= (x[n - nmin] * tmpy[n - nmin] * (pxca[n - nmin] - pysa[n - nmin]) +
0284                            y[n - nmin] * (1.f + tmpy[n - nmin] * (pyca[n - nmin] + pxsa[n - nmin]))) *
0285                           oor0[n - nmin];
0286       }
0287 #pragma omp simd
0288       for (int n = nmin; n < nmax; ++n) {
0289         //now r0 depends on ipt and phi as well
0290         tmp[n - nmin] = dadipt[n - nmin] * ipt[n - nmin];
0291       }
0292 #pragma omp simd
0293       for (int n = nmin; n < nmax; ++n) {
0294         dDdipt[n - nmin] -= k[n - nmin] *
0295                             (x[n - nmin] * (pxca[n - nmin] * tmp[n - nmin] - pysa[n - nmin] * tmp[n - nmin] -
0296                                             pyca[n - nmin] - pxsa[n - nmin] + pyin[n - nmin]) +
0297                              y[n - nmin] * (pyca[n - nmin] * tmp[n - nmin] + pxsa[n - nmin] * tmp[n - nmin] -
0298                                             pysa[n - nmin] + pxca[n - nmin] - pxin[n - nmin])) *
0299                             pt[n - nmin] * oor0[n - nmin];
0300       }
0301 #pragma omp simd
0302       for (int n = nmin; n < nmax; ++n) {
0303         dDdphi[n - nmin] += k[n - nmin] *
0304                             (x[n - nmin] * (pysa[n - nmin] - pxin[n - nmin] + pxca[n - nmin]) -
0305                              y[n - nmin] * (pxsa[n - nmin] - pyin[n - nmin] + pyca[n - nmin])) *
0306                             oor0[n - nmin];
0307       }
0308     }
0309 
0310 #pragma omp simd
0311     for (int n = nmin; n < nmax; ++n) {
0312       //update parameters
0313       outPar(n, 0, 0) = outPar(n, 0, 0) + 2.f * k[n - nmin] * sinah[n - nmin] *
0314                                               (pxin[n - nmin] * cosah[n - nmin] - pyin[n - nmin] * sinah[n - nmin]);
0315       outPar(n, 1, 0) = outPar(n, 1, 0) + 2.f * k[n - nmin] * sinah[n - nmin] *
0316                                               (pyin[n - nmin] * cosah[n - nmin] + pxin[n - nmin] * sinah[n - nmin]);
0317       pxinold[n - nmin] = pxin[n - nmin];  //copy before overwriting
0318       pxin[n - nmin] = pxin[n - nmin] * cosa[n - nmin] - pyin[n - nmin] * sina[n - nmin];
0319       pyin[n - nmin] = pyin[n - nmin] * cosa[n - nmin] + pxinold[n - nmin] * sina[n - nmin];
0320     }
0321     for (int n = nmin; n < nmax; ++n) {
0322       dprint_np(n,
0323                 "outPar(n, 0, 0)=" << outPar(n, 0, 0) << " outPar(n, 1, 0)=" << outPar(n, 1, 0)
0324                                    << " pxin=" << pxin[n - nmin] << " pyin=" << pyin[n - nmin]);
0325     }
0326   }  // iteration loop
0327 
0328   float alpha[nmax - nmin];
0329   float dadphi[nmax - nmin];
0330 
0331 #pragma omp simd
0332   for (int n = nmin; n < nmax; ++n) {
0333     alpha[n - nmin] = D[n - nmin] * ipt[n - nmin] * kinv[n - nmin];
0334     dadx[n - nmin] = dDdx[n - nmin] * ipt[n - nmin] * kinv[n - nmin];
0335     dady[n - nmin] = dDdy[n - nmin] * ipt[n - nmin] * kinv[n - nmin];
0336     dadipt[n - nmin] = (ipt[n - nmin] * dDdipt[n - nmin] + D[n - nmin]) * kinv[n - nmin];
0337     dadphi[n - nmin] = dDdphi[n - nmin] * ipt[n - nmin] * kinv[n - nmin];
0338   }
0339 
0340   if constexpr (Config::useTrigApprox) {
0341 #pragma omp simd
0342     for (int n = nmin; n < nmax; ++n) {
0343       sincos4(alpha[n - nmin], sina[n - nmin], cosa[n - nmin]);
0344     }
0345   } else {
0346 #pragma omp simd
0347     for (int n = nmin; n < nmax; ++n) {
0348       cosa[n - nmin] = std::cos(alpha[n - nmin]);
0349       sina[n - nmin] = std::sin(alpha[n - nmin]);
0350     }
0351   }
0352 #pragma omp simd
0353   for (int n = nmin; n < nmax; ++n) {
0354     errorProp(n, 0, 0) = 1.f + k[n - nmin] * dadx[n - nmin] *
0355                                    (cosPorT[n - nmin] * cosa[n - nmin] - sinPorT[n - nmin] * sina[n - nmin]) *
0356                                    pt[n - nmin];
0357     errorProp(n, 0, 1) = k[n - nmin] * dady[n - nmin] *
0358                          (cosPorT[n - nmin] * cosa[n - nmin] - sinPorT[n - nmin] * sina[n - nmin]) * pt[n - nmin];
0359     errorProp(n, 0, 2) = 0.f;
0360     errorProp(n, 0, 3) =
0361         k[n - nmin] *
0362         (cosPorT[n - nmin] * (ipt[n - nmin] * dadipt[n - nmin] * cosa[n - nmin] - sina[n - nmin]) +
0363          sinPorT[n - nmin] * ((1.f - cosa[n - nmin]) - ipt[n - nmin] * dadipt[n - nmin] * sina[n - nmin])) *
0364         pt[n - nmin] * pt[n - nmin];
0365     errorProp(n, 0, 4) =
0366         k[n - nmin] *
0367         (cosPorT[n - nmin] * dadphi[n - nmin] * cosa[n - nmin] - sinPorT[n - nmin] * dadphi[n - nmin] * sina[n - nmin] -
0368          sinPorT[n - nmin] * sina[n - nmin] + cosPorT[n - nmin] * cosa[n - nmin] - cosPorT[n - nmin]) *
0369         pt[n - nmin];
0370     errorProp(n, 0, 5) = 0.f;
0371   }
0372 
0373 #pragma omp simd
0374   for (int n = nmin; n < nmax; ++n) {
0375     errorProp(n, 1, 0) = k[n - nmin] * dadx[n - nmin] *
0376                          (sinPorT[n - nmin] * cosa[n - nmin] + cosPorT[n - nmin] * sina[n - nmin]) * pt[n - nmin];
0377     errorProp(n, 1, 1) = 1.f + k[n - nmin] * dady[n - nmin] *
0378                                    (sinPorT[n - nmin] * cosa[n - nmin] + cosPorT[n - nmin] * sina[n - nmin]) *
0379                                    pt[n - nmin];
0380     errorProp(n, 1, 2) = 0.f;
0381     errorProp(n, 1, 3) =
0382         k[n - nmin] *
0383         (sinPorT[n - nmin] * (ipt[n - nmin] * dadipt[n - nmin] * cosa[n - nmin] - sina[n - nmin]) +
0384          cosPorT[n - nmin] * (ipt[n - nmin] * dadipt[n - nmin] * sina[n - nmin] - (1.f - cosa[n - nmin]))) *
0385         pt[n - nmin] * pt[n - nmin];
0386     errorProp(n, 1, 4) =
0387         k[n - nmin] *
0388         (sinPorT[n - nmin] * dadphi[n - nmin] * cosa[n - nmin] + cosPorT[n - nmin] * dadphi[n - nmin] * sina[n - nmin] +
0389          sinPorT[n - nmin] * cosa[n - nmin] + cosPorT[n - nmin] * sina[n - nmin] - sinPorT[n - nmin]) *
0390         pt[n - nmin];
0391     errorProp(n, 1, 5) = 0.f;
0392   }
0393 
0394 #pragma omp simd
0395   for (int n = nmin; n < nmax; ++n) {
0396     //no trig approx here, theta can be large
0397     cosPorT[n - nmin] = std::cos(theta[n - nmin]);
0398   }
0399 
0400 #pragma omp simd
0401   for (int n = nmin; n < nmax; ++n) {
0402     sinPorT[n - nmin] = std::sin(theta[n - nmin]);
0403   }
0404 
0405 #pragma omp simd
0406   for (int n = nmin; n < nmax; ++n) {
0407     //redefine sinPorT as 1./sinPorT to reduce the number of temporaries
0408     sinPorT[n - nmin] = 1.f / sinPorT[n - nmin];
0409   }
0410 
0411 #pragma omp simd
0412   for (int n = nmin; n < nmax; ++n) {
0413     outPar(n, 2, 0) =
0414         inPar(n, 2, 0) + k[n - nmin] * alpha[n - nmin] * cosPorT[n - nmin] * pt[n - nmin] * sinPorT[n - nmin];
0415     errorProp(n, 2, 0) = k[n - nmin] * cosPorT[n - nmin] * dadx[n - nmin] * pt[n - nmin] * sinPorT[n - nmin];
0416     errorProp(n, 2, 1) = k[n - nmin] * cosPorT[n - nmin] * dady[n - nmin] * pt[n - nmin] * sinPorT[n - nmin];
0417     errorProp(n, 2, 2) = 1.f;
0418     errorProp(n, 2, 3) = k[n - nmin] * cosPorT[n - nmin] * (ipt[n - nmin] * dadipt[n - nmin] - alpha[n - nmin]) *
0419                          pt[n - nmin] * pt[n - nmin] * sinPorT[n - nmin];
0420     errorProp(n, 2, 4) = k[n - nmin] * dadphi[n - nmin] * cosPorT[n - nmin] * pt[n - nmin] * sinPorT[n - nmin];
0421     errorProp(n, 2, 5) = -k[n - nmin] * alpha[n - nmin] * pt[n - nmin] * sinPorT[n - nmin] * sinPorT[n - nmin];
0422   }
0423 
0424 #pragma omp simd
0425   for (int n = nmin; n < nmax; ++n) {
0426     outPar(n, 3, 0) = ipt[n - nmin];
0427     errorProp(n, 3, 0) = 0.f;
0428     errorProp(n, 3, 1) = 0.f;
0429     errorProp(n, 3, 2) = 0.f;
0430     errorProp(n, 3, 3) = 1.f;
0431     errorProp(n, 3, 4) = 0.f;
0432     errorProp(n, 3, 5) = 0.f;
0433   }
0434 
0435 #pragma omp simd
0436   for (int n = nmin; n < nmax; ++n) {
0437     outPar(n, 4, 0) = inPar(n, 4, 0) + alpha[n - nmin];
0438     errorProp(n, 4, 0) = dadx[n - nmin];
0439     errorProp(n, 4, 1) = dady[n - nmin];
0440     errorProp(n, 4, 2) = 0.f;
0441     errorProp(n, 4, 3) = dadipt[n - nmin];
0442     errorProp(n, 4, 4) = 1.f + dadphi[n - nmin];
0443     errorProp(n, 4, 5) = 0.f;
0444   }
0445 
0446 #pragma omp simd
0447   for (int n = nmin; n < nmax; ++n) {
0448     outPar(n, 5, 0) = theta[n - nmin];
0449     errorProp(n, 5, 0) = 0.f;
0450     errorProp(n, 5, 1) = 0.f;
0451     errorProp(n, 5, 2) = 0.f;
0452     errorProp(n, 5, 3) = 0.f;
0453     errorProp(n, 5, 4) = 0.f;
0454     errorProp(n, 5, 5) = 1.f;
0455   }
0456 
0457   for (int n = nmin; n < nmax; ++n) {
0458     dprint_np(n,
0459               "propagation end, dump parameters\n"
0460                   << "   pos = " << outPar(n, 0, 0) << " " << outPar(n, 1, 0) << " " << outPar(n, 2, 0) << "\t\t r="
0461                   << std::sqrt(outPar(n, 0, 0) * outPar(n, 0, 0) + outPar(n, 1, 0) * outPar(n, 1, 0)) << std::endl
0462                   << "   mom = " << std::cos(outPar(n, 4, 0)) / outPar(n, 3, 0) << " "
0463                   << std::sin(outPar(n, 4, 0)) / outPar(n, 3, 0) << " " << 1. / (outPar(n, 3, 0) * tan(outPar(n, 5, 0)))
0464                   << "\t\tpT=" << 1. / std::abs(outPar(n, 3, 0)) << std::endl);
0465   }
0466 
0467 #ifdef DEBUG
0468   for (int n = nmin; n < nmax; ++n) {
0469     if (debug && g_debug && n < N_proc) {
0470       dmutex_guard;
0471       std::cout << n << ": jacobian" << std::endl;
0472       printf("%5f %5f %5f %5f %5f %5f\n",
0473              errorProp(n, 0, 0),
0474              errorProp(n, 0, 1),
0475              errorProp(n, 0, 2),
0476              errorProp(n, 0, 3),
0477              errorProp(n, 0, 4),
0478              errorProp(n, 0, 5));
0479       printf("%5f %5f %5f %5f %5f %5f\n",
0480              errorProp(n, 1, 0),
0481              errorProp(n, 1, 1),
0482              errorProp(n, 1, 2),
0483              errorProp(n, 1, 3),
0484              errorProp(n, 1, 4),
0485              errorProp(n, 1, 5));
0486       printf("%5f %5f %5f %5f %5f %5f\n",
0487              errorProp(n, 2, 0),
0488              errorProp(n, 2, 1),
0489              errorProp(n, 2, 2),
0490              errorProp(n, 2, 3),
0491              errorProp(n, 2, 4),
0492              errorProp(n, 2, 5));
0493       printf("%5f %5f %5f %5f %5f %5f\n",
0494              errorProp(n, 3, 0),
0495              errorProp(n, 3, 1),
0496              errorProp(n, 3, 2),
0497              errorProp(n, 3, 3),
0498              errorProp(n, 3, 4),
0499              errorProp(n, 3, 5));
0500       printf("%5f %5f %5f %5f %5f %5f\n",
0501              errorProp(n, 4, 0),
0502              errorProp(n, 4, 1),
0503              errorProp(n, 4, 2),
0504              errorProp(n, 4, 3),
0505              errorProp(n, 4, 4),
0506              errorProp(n, 4, 5));
0507       printf("%5f %5f %5f %5f %5f %5f\n",
0508              errorProp(n, 5, 0),
0509              errorProp(n, 5, 1),
0510              errorProp(n, 5, 2),
0511              errorProp(n, 5, 3),
0512              errorProp(n, 5, 4),
0513              errorProp(n, 5, 5));
0514       printf("\n");
0515     }
0516   }
0517 #endif
0518 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!