CMS 3D CMS Logo
PropagationMPlex.cc
Go to the documentation of this file.
1 #include "FWCore/Utilities/interface/CMSUnrollLoop.h"
2 #include "RecoTracker/MkFitCore/interface/PropagationConfig.h"
3 #include "RecoTracker/MkFitCore/interface/Config.h"
4 #include "RecoTracker/MkFitCore/interface/TrackerInfo.h"
5 
6 #include "PropagationMPlex.h"
7 
8 //#define DEBUG
9 #include "Debug.h"
10 
11 namespace mkfit {
12 
13  //==============================================================================
14  // propagateLineToRMPlex
15  //==============================================================================
16 
17  void propagateLineToRMPlex(const MPlexLS& psErr,
18  const MPlexLV& psPar,
19  const MPlexHS& msErr,
20  const MPlexHV& msPar,
21  MPlexLS& outErr,
22  MPlexLV& outPar,
23  const int N_proc) {
24  // XXX Regenerate parts below with a script.
25 
26  const Matriplex::idx_t N = NN;
27 
28 #pragma omp simd
29  for (int n = 0; n < NN; ++n) {
30  const float cosA = (psPar[0 * N + n] * psPar[3 * N + n] + psPar[1 * N + n] * psPar[4 * N + n]) /
31  (std::sqrt((psPar[0 * N + n] * psPar[0 * N + n] + psPar[1 * N + n] * psPar[1 * N + n]) *
32  (psPar[3 * N + n] * psPar[3 * N + n] + psPar[4 * N + n] * psPar[4 * N + n])));
33  const float dr = (hipo(msPar[0 * N + n], msPar[1 * N + n]) - hipo(psPar[0 * N + n], psPar[1 * N + n])) / cosA;
34 
35  dprint_np(n, "propagateLineToRMPlex dr=" << dr);
36 
37  const float pt = hipo(psPar[3 * N + n], psPar[4 * N + n]);
38  const float p = dr / pt; // path
39  const float psq = p * p;
40 
41  outPar[0 * N + n] = psPar[0 * N + n] + p * psPar[3 * N + n];
42  outPar[1 * N + n] = psPar[1 * N + n] + p * psPar[4 * N + n];
43  outPar[2 * N + n] = psPar[2 * N + n] + p * psPar[5 * N + n];
44  outPar[3 * N + n] = psPar[3 * N + n];
45  outPar[4 * N + n] = psPar[4 * N + n];
46  outPar[5 * N + n] = psPar[5 * N + n];
47 
48  {
49  const MPlexLS& A = psErr;
50  MPlexLS& B = outErr;
51 
52  B.fArray[0 * N + n] = A.fArray[0 * N + n];
53  B.fArray[1 * N + n] = A.fArray[1 * N + n];
54  B.fArray[2 * N + n] = A.fArray[2 * N + n];
55  B.fArray[3 * N + n] = A.fArray[3 * N + n];
56  B.fArray[4 * N + n] = A.fArray[4 * N + n];
57  B.fArray[5 * N + n] = A.fArray[5 * N + n];
58  B.fArray[6 * N + n] = A.fArray[6 * N + n] + p * A.fArray[0 * N + n];
59  B.fArray[7 * N + n] = A.fArray[7 * N + n] + p * A.fArray[1 * N + n];
60  B.fArray[8 * N + n] = A.fArray[8 * N + n] + p * A.fArray[3 * N + n];
61  B.fArray[9 * N + n] =
62  A.fArray[9 * N + n] + p * (A.fArray[6 * N + n] + A.fArray[6 * N + n]) + psq * A.fArray[0 * N + n];
63  B.fArray[10 * N + n] = A.fArray[10 * N + n] + p * A.fArray[1 * N + n];
64  B.fArray[11 * N + n] = A.fArray[11 * N + n] + p * A.fArray[2 * N + n];
65  B.fArray[12 * N + n] = A.fArray[12 * N + n] + p * A.fArray[4 * N + n];
66  B.fArray[13 * N + n] =
67  A.fArray[13 * N + n] + p * (A.fArray[7 * N + n] + A.fArray[10 * N + n]) + psq * A.fArray[1 * N + n];
68  B.fArray[14 * N + n] =
69  A.fArray[14 * N + n] + p * (A.fArray[11 * N + n] + A.fArray[11 * N + n]) + psq * A.fArray[2 * N + n];
70  B.fArray[15 * N + n] = A.fArray[15 * N + n] + p * A.fArray[3 * N + n];
71  B.fArray[16 * N + n] = A.fArray[16 * N + n] + p * A.fArray[4 * N + n];
72  B.fArray[17 * N + n] = A.fArray[17 * N + n] + p * A.fArray[5 * N + n];
73  B.fArray[18 * N + n] =
74  A.fArray[18 * N + n] + p * (A.fArray[8 * N + n] + A.fArray[15 * N + n]) + psq * A.fArray[3 * N + n];
75  B.fArray[19 * N + n] =
76  A.fArray[19 * N + n] + p * (A.fArray[12 * N + n] + A.fArray[16 * N + n]) + psq * A.fArray[4 * N + n];
77  B.fArray[20 * N + n] =
78  A.fArray[20 * N + n] + p * (A.fArray[17 * N + n] + A.fArray[17 * N + n]) + psq * A.fArray[5 * N + n];
79  }
80 
81  dprint_np(n, "propagateLineToRMPlex arrive at r=" << hipo(outPar[0 * N + n], outPar[1 * N + n]));
82  }
83  }
84 
85 } // end namespace mkfit
86 
87 //==============================================================================
88 // propagateHelixToRMPlex
89 //==============================================================================
90 
91 namespace {
92  using namespace mkfit;
93 
94  void MultHelixProp(const MPlexLL& A, const MPlexLS& B, MPlexLL& C) {
95  // C = A * B
96 
97  typedef float T;
98  const Matriplex::idx_t N = NN;
99 
100  const T* a = A.fArray;
101  ASSUME_ALIGNED(a, 64);
102  const T* b = B.fArray;
103  ASSUME_ALIGNED(b, 64);
104  T* c = C.fArray;
105  ASSUME_ALIGNED(c, 64);
106 
107 #include "MultHelixProp.ah"
108  }
109 
110  void MultHelixPropTransp(const MPlexLL& A, const MPlexLL& B, MPlexLS& C) {
111  // C = B * AT;
112 
113  typedef float T;
114  const Matriplex::idx_t N = NN;
115 
116  const T* a = A.fArray;
117  ASSUME_ALIGNED(a, 64);
118  const T* b = B.fArray;
119  ASSUME_ALIGNED(b, 64);
120  T* c = C.fArray;
121  ASSUME_ALIGNED(c, 64);
122 
123 #include "MultHelixPropTransp.ah"
124  }
125 
126  void MultHelixPropTemp(const MPlexLL& A, const MPlexLL& B, MPlexLL& C, int n) {
127  // C = A * B
128 
129  typedef float T;
130  const Matriplex::idx_t N = NN;
131 
132  const T* a = A.fArray;
133  ASSUME_ALIGNED(a, 64);
134  const T* b = B.fArray;
135  ASSUME_ALIGNED(b, 64);
136  T* c = C.fArray;
137  ASSUME_ALIGNED(c, 64);
138 
139  c[0 * N + n] = a[0 * N + n] * b[0 * N + n] + a[1 * N + n] * b[6 * N + n] + a[2 * N + n] * b[12 * N + n] +
140  a[4 * N + n] * b[24 * N + n];
141  c[1 * N + n] = a[0 * N + n] * b[1 * N + n] + a[1 * N + n] * b[7 * N + n] + a[2 * N + n] * b[13 * N + n] +
142  a[4 * N + n] * b[25 * N + n];
143  c[2 * N + n] = a[2 * N + n];
144  c[3 * N + n] = a[0 * N + n] * b[3 * N + n] + a[1 * N + n] * b[9 * N + n] + a[2 * N + n] * b[15 * N + n] +
145  a[3 * N + n] + a[4 * N + n] * b[27 * N + n];
146  c[4 * N + n] = a[0 * N + n] * b[4 * N + n] + a[1 * N + n] * b[10 * N + n] + a[4 * N + n];
147  c[5 * N + n] = a[2 * N + n] * b[17 * N + n] + a[5 * N + n];
148  c[6 * N + n] = a[6 * N + n] * b[0 * N + n] + a[7 * N + n] * b[6 * N + n] + a[8 * N + n] * b[12 * N + n] +
149  a[10 * N + n] * b[24 * N + n];
150  c[7 * N + n] = a[6 * N + n] * b[1 * N + n] + a[7 * N + n] * b[7 * N + n] + a[8 * N + n] * b[13 * N + n] +
151  a[10 * N + n] * b[25 * N + n];
152  c[8 * N + n] = a[8 * N + n];
153  c[9 * N + n] = a[6 * N + n] * b[3 * N + n] + a[7 * N + n] * b[9 * N + n] + a[8 * N + n] * b[15 * N + n] +
154  a[9 * N + n] + a[10 * N + n] * b[27 * N + n];
155  c[10 * N + n] = a[6 * N + n] * b[4 * N + n] + a[7 * N + n] * b[10 * N + n] + a[10 * N + n];
156  c[11 * N + n] = a[8 * N + n] * b[17 * N + n] + a[11 * N + n];
157  c[12 * N + n] = a[12 * N + n] * b[0 * N + n] + a[13 * N + n] * b[6 * N + n] + a[14 * N + n] * b[12 * N + n] +
158  a[16 * N + n] * b[24 * N + n];
159  c[13 * N + n] = a[12 * N + n] * b[1 * N + n] + a[13 * N + n] * b[7 * N + n] + a[14 * N + n] * b[13 * N + n] +
160  a[16 * N + n] * b[25 * N + n];
161  c[14 * N + n] = a[14 * N + n];
162  c[15 * N + n] = a[12 * N + n] * b[3 * N + n] + a[13 * N + n] * b[9 * N + n] + a[14 * N + n] * b[15 * N + n] +
163  a[15 * N + n] + a[16 * N + n] * b[27 * N + n];
164  c[16 * N + n] = a[12 * N + n] * b[4 * N + n] + a[13 * N + n] * b[10 * N + n] + a[16 * N + n];
165  c[17 * N + n] = a[14 * N + n] * b[17 * N + n] + a[17 * N + n];
166  c[18 * N + n] = a[18 * N + n] * b[0 * N + n] + a[19 * N + n] * b[6 * N + n] + a[20 * N + n] * b[12 * N + n] +
167  a[22 * N + n] * b[24 * N + n];
168  c[19 * N + n] = a[18 * N + n] * b[1 * N + n] + a[19 * N + n] * b[7 * N + n] + a[20 * N + n] * b[13 * N + n] +
169  a[22 * N + n] * b[25 * N + n];
170  c[20 * N + n] = a[20 * N + n];
171  c[21 * N + n] = a[18 * N + n] * b[3 * N + n] + a[19 * N + n] * b[9 * N + n] + a[20 * N + n] * b[15 * N + n] +
172  a[21 * N + n] + a[22 * N + n] * b[27 * N + n];
173  c[22 * N + n] = a[18 * N + n] * b[4 * N + n] + a[19 * N + n] * b[10 * N + n] + a[22 * N + n];
174  c[23 * N + n] = a[20 * N + n] * b[17 * N + n] + a[23 * N + n];
175  c[24 * N + n] = a[24 * N + n] * b[0 * N + n] + a[25 * N + n] * b[6 * N + n] + a[26 * N + n] * b[12 * N + n] +
176  a[28 * N + n] * b[24 * N + n];
177  c[25 * N + n] = a[24 * N + n] * b[1 * N + n] + a[25 * N + n] * b[7 * N + n] + a[26 * N + n] * b[13 * N + n] +
178  a[28 * N + n] * b[25 * N + n];
179  c[26 * N + n] = a[26 * N + n];
180  c[27 * N + n] = a[24 * N + n] * b[3 * N + n] + a[25 * N + n] * b[9 * N + n] + a[26 * N + n] * b[15 * N + n] +
181  a[27 * N + n] + a[28 * N + n] * b[27 * N + n];
182  c[28 * N + n] = a[24 * N + n] * b[4 * N + n] + a[25 * N + n] * b[10 * N + n] + a[28 * N + n];
183  c[29 * N + n] = a[26 * N + n] * b[17 * N + n] + a[29 * N + n];
184  c[30 * N + n] = a[30 * N + n] * b[0 * N + n] + a[31 * N + n] * b[6 * N + n] + a[32 * N + n] * b[12 * N + n] +
185  a[34 * N + n] * b[24 * N + n];
186  c[31 * N + n] = a[30 * N + n] * b[1 * N + n] + a[31 * N + n] * b[7 * N + n] + a[32 * N + n] * b[13 * N + n] +
187  a[34 * N + n] * b[25 * N + n];
188  c[32 * N + n] = a[32 * N + n];
189  c[33 * N + n] = a[30 * N + n] * b[3 * N + n] + a[31 * N + n] * b[9 * N + n] + a[32 * N + n] * b[15 * N + n] +
190  a[33 * N + n] + a[34 * N + n] * b[27 * N + n];
191  c[34 * N + n] = a[30 * N + n] * b[4 * N + n] + a[31 * N + n] * b[10 * N + n] + a[34 * N + n];
192  c[35 * N + n] = a[32 * N + n] * b[17 * N + n] + a[35 * N + n];
193  }
194 
195 } // end unnamed namespace
196 
197 //==============================================================================
198 
199 namespace mkfit {
200 
201  void helixAtRFromIterativeCCSFullJac(const MPlexLV& inPar,
202  const MPlexQI& inChg,
203  const MPlexQF& msRad,
204  MPlexLV& outPar,
205  MPlexLL& errorProp,
206  const int N_proc) {
207  errorProp.setVal(0.f);
208  MPlexLL errorPropTmp(0.f); //initialize to zero
209  MPlexLL errorPropSwap(0.f); //initialize to zero
210 
211  // loop does not vectorize with llvm16, and it issues a warning
212  // that apparently can't be suppressed with a pragma. Needs to
213  // be rechecked if future llvm versions improve vectorization.
214 #if !defined(__clang__)
215 #pragma omp simd
216 #endif
217  for (int n = 0; n < NN; ++n) {
218  //initialize erroProp to identity matrix
219  errorProp(n, 0, 0) = 1.f;
220  errorProp(n, 1, 1) = 1.f;
221  errorProp(n, 2, 2) = 1.f;
222  errorProp(n, 3, 3) = 1.f;
223  errorProp(n, 4, 4) = 1.f;
224  errorProp(n, 5, 5) = 1.f;
225 
226  const float k = inChg.constAt(n, 0, 0) * 100.f / (-Const::sol * Config::Bfield);
227  const float r = msRad.constAt(n, 0, 0);
228  float r0 = hipo(inPar.constAt(n, 0, 0), inPar.constAt(n, 1, 0));
229 
230  if (std::abs(r - r0) < 0.0001f) {
231  dprint_np(n, "distance less than 1mum, skip");
232  continue;
233  }
234 
235  const float ipt = inPar.constAt(n, 3, 0);
236  const float phiin = inPar.constAt(n, 4, 0);
237  const float theta = inPar.constAt(n, 5, 0);
238 
239  //set those that are 1. before iterations
240  errorPropTmp(n, 2, 2) = 1.f;
241  errorPropTmp(n, 3, 3) = 1.f;
242  errorPropTmp(n, 4, 4) = 1.f;
243  errorPropTmp(n, 5, 5) = 1.f;
244 
245  float cosah = 0., sinah = 0.;
246  //no trig approx here, phi and theta can be large
247  float cosP = std::cos(phiin), sinP = std::sin(phiin);
248  const float cosT = std::cos(theta), sinT = std::sin(theta);
249  float pxin = cosP / ipt;
250  float pyin = sinP / ipt;
251 
252  CMS_UNROLL_LOOP_COUNT(Config::Niter)
253  for (int i = 0; i < Config::Niter; ++i) {
254  dprint_np(n,
255  std::endl
256  << "attempt propagation from r=" << r0 << " to r=" << r << std::endl
257  << "x=" << outPar.At(n, 0, 0) << " y=" << outPar.At(n, 1, 0) << " z=" << outPar.At(n, 2, 0)
258  << " px=" << std::cos(phiin) / ipt << " py=" << std::sin(phiin) / ipt
259  << " pz=" << 1.f / (ipt * tan(theta)) << " q=" << inChg.constAt(n, 0, 0) << std::endl);
260 
261  r0 = hipo(outPar.constAt(n, 0, 0), outPar.constAt(n, 1, 0));
262  const float ialpha = (r - r0) * ipt / k;
263  //alpha+=ialpha;
264 
265  if constexpr (Config::useTrigApprox) {
266  sincos4(ialpha * 0.5f, sinah, cosah);
267  } else {
268  cosah = std::cos(ialpha * 0.5f);
269  sinah = std::sin(ialpha * 0.5f);
270  }
271  const float cosa = 1.f - 2.f * sinah * sinah;
272  const float sina = 2.f * sinah * cosah;
273 
274  //derivatives of alpha
275  const float dadx = -outPar.At(n, 0, 0) * ipt / (k * r0);
276  const float dady = -outPar.At(n, 1, 0) * ipt / (k * r0);
277  const float dadipt = (r - r0) / k;
278 
279  outPar.At(n, 0, 0) = outPar.constAt(n, 0, 0) + 2.f * k * sinah * (pxin * cosah - pyin * sinah);
280  outPar.At(n, 1, 0) = outPar.constAt(n, 1, 0) + 2.f * k * sinah * (pyin * cosah + pxin * sinah);
281  const float pxinold = pxin; //copy before overwriting
282  pxin = pxin * cosa - pyin * sina;
283  pyin = pyin * cosa + pxinold * sina;
284 
285  //need phi at origin, so this goes before redefining phi
286  //no trig approx here, phi can be large
287  cosP = std::cos(outPar.At(n, 4, 0));
288  sinP = std::sin(outPar.At(n, 4, 0));
289 
290  outPar.At(n, 2, 0) = outPar.constAt(n, 2, 0) + k * ialpha * cosT / (ipt * sinT);
291  outPar.At(n, 3, 0) = ipt;
292  outPar.At(n, 4, 0) = outPar.constAt(n, 4, 0) + ialpha;
293  outPar.At(n, 5, 0) = theta;
294 
295  errorPropTmp(n, 0, 0) = 1.f + k * (cosP * dadx * cosa - sinP * dadx * sina) / ipt;
296  errorPropTmp(n, 0, 1) = k * (cosP * dady * cosa - sinP * dady * sina) / ipt;
297  errorPropTmp(n, 0, 3) =
298  k * (cosP * (ipt * dadipt * cosa - sina) + sinP * ((1.f - cosa) - ipt * dadipt * sina)) / (ipt * ipt);
299  errorPropTmp(n, 0, 4) = -k * (sinP * sina + cosP * (1.f - cosa)) / ipt;
300 
301  errorPropTmp(n, 1, 0) = k * (sinP * dadx * cosa + cosP * dadx * sina) / ipt;
302  errorPropTmp(n, 1, 1) = 1.f + k * (sinP * dady * cosa + cosP * dady * sina) / ipt;
303  errorPropTmp(n, 1, 3) =
304  k * (sinP * (ipt * dadipt * cosa - sina) + cosP * (ipt * dadipt * sina - (1.f - cosa))) / (ipt * ipt);
305  errorPropTmp(n, 1, 4) = k * (cosP * sina - sinP * (1.f - cosa)) / ipt;
306 
307  errorPropTmp(n, 2, 0) = k * cosT * dadx / (ipt * sinT);
308  errorPropTmp(n, 2, 1) = k * cosT * dady / (ipt * sinT);
309  errorPropTmp(n, 2, 3) = k * cosT * (ipt * dadipt - ialpha) / (ipt * ipt * sinT);
310  errorPropTmp(n, 2, 5) = -k * ialpha / (ipt * sinT * sinT);
311 
312  errorPropTmp(n, 4, 0) = dadx;
313  errorPropTmp(n, 4, 1) = dady;
314  errorPropTmp(n, 4, 3) = dadipt;
315 
316  MultHelixPropTemp(errorProp, errorPropTmp, errorPropSwap, n);
317  errorProp = errorPropSwap;
318  }
319 
320  dprint_np(
321  n,
322  "propagation end, dump parameters"
323  << std::endl
324  << "pos = " << outPar.At(n, 0, 0) << " " << outPar.At(n, 1, 0) << " " << outPar.At(n, 2, 0) << std::endl
325  << "mom = " << std::cos(outPar.At(n, 4, 0)) / outPar.At(n, 3, 0) << " "
326  << std::sin(outPar.At(n, 4, 0)) / outPar.At(n, 3, 0) << " "
327  << 1. / (outPar.At(n, 3, 0) * tan(outPar.At(n, 5, 0)))
328  << " r=" << std::sqrt(outPar.At(n, 0, 0) * outPar.At(n, 0, 0) + outPar.At(n, 1, 0) * outPar.At(n, 1, 0))
329  << " pT=" << 1. / std::abs(outPar.At(n, 3, 0)) << std::endl);
330 
331 #ifdef DEBUG
332  if (debug && g_debug && n < N_proc) {
333  dmutex_guard;
334  std::cout << n << " jacobian" << std::endl;
335  printf("%5f %5f %5f %5f %5f %5f\n",
336  errorProp(n, 0, 0),
337  errorProp(n, 0, 1),
338  errorProp(n, 0, 2),
339  errorProp(n, 0, 3),
340  errorProp(n, 0, 4),
341  errorProp(n, 0, 5));
342  printf("%5f %5f %5f %5f %5f %5f\n",
343  errorProp(n, 1, 0),
344  errorProp(n, 1, 1),
345  errorProp(n, 1, 2),
346  errorProp(n, 1, 3),
347  errorProp(n, 1, 4),
348  errorProp(n, 1, 5));
349  printf("%5f %5f %5f %5f %5f %5f\n",
350  errorProp(n, 2, 0),
351  errorProp(n, 2, 1),
352  errorProp(n, 2, 2),
353  errorProp(n, 2, 3),
354  errorProp(n, 2, 4),
355  errorProp(n, 2, 5));
356  printf("%5f %5f %5f %5f %5f %5f\n",
357  errorProp(n, 3, 0),
358  errorProp(n, 3, 1),
359  errorProp(n, 3, 2),
360  errorProp(n, 3, 3),
361  errorProp(n, 3, 4),
362  errorProp(n, 3, 5));
363  printf("%5f %5f %5f %5f %5f %5f\n",
364  errorProp(n, 4, 0),
365  errorProp(n, 4, 1),
366  errorProp(n, 4, 2),
367  errorProp(n, 4, 3),
368  errorProp(n, 4, 4),
369  errorProp(n, 4, 5));
370  printf("%5f %5f %5f %5f %5f %5f\n",
371  errorProp(n, 5, 0),
372  errorProp(n, 5, 1),
373  errorProp(n, 5, 2),
374  errorProp(n, 5, 3),
375  errorProp(n, 5, 4),
376  errorProp(n, 5, 5));
377  }
378 #endif
379  }
380  }
381 
382 } // end namespace mkfit
383 
384 // ============================================================================
385 // BEGIN STUFF FROM PropagationMPlex.icc
386 
387 namespace {
388 
389  //========================================================================================
390  // helixAtR
391  //========================================================================================
392 
393  void helixAtRFromIterativeCCS_impl(const MPlexLV& __restrict__ inPar,
394  const MPlexQI& __restrict__ inChg,
395  const MPlexQF& __restrict__ msRad,
396  MPlexLV& __restrict__ outPar,
397  MPlexLL& __restrict__ errorProp,
398  MPlexQI& __restrict__ outFailFlag, // expected to be initialized to 0
399  const int nmin,
400  const int nmax,
401  const int N_proc,
402  const PropagationFlags& pf) {
403  // bool debug = true;
404 
405 #pragma omp simd
406  for (int n = nmin; n < nmax; ++n) {
407  //initialize erroProp to identity matrix
408  errorProp(n, 0, 0) = 1.f;
409  errorProp(n, 1, 1) = 1.f;
410  errorProp(n, 2, 2) = 1.f;
411  errorProp(n, 3, 3) = 1.f;
412  errorProp(n, 4, 4) = 1.f;
413  errorProp(n, 5, 5) = 1.f;
414  }
415  float r0[nmax - nmin];
416 #pragma omp simd
417  for (int n = nmin; n < nmax; ++n) {
418  //initialize erroProp to identity matrix
419  r0[n - nmin] = hipo(inPar(n, 0, 0), inPar(n, 1, 0));
420  }
421  float k[nmax - nmin];
422  if (pf.use_param_b_field) {
423 #pragma omp simd
424  for (int n = nmin; n < nmax; ++n) {
425  k[n - nmin] = inChg(n, 0, 0) * 100.f / (-Const::sol * Config::bFieldFromZR(inPar(n, 2, 0), r0[n - nmin]));
426  }
427  } else {
428 #pragma omp simd
429  for (int n = nmin; n < nmax; ++n) {
430  k[n - nmin] = inChg(n, 0, 0) * 100.f / (-Const::sol * Config::Bfield);
431  }
432  }
433  float r[nmax - nmin];
434 #pragma omp simd
435  for (int n = nmin; n < nmax; ++n) {
436  r[n - nmin] = msRad(n, 0, 0);
437  }
438  float xin[nmax - nmin];
439  float yin[nmax - nmin];
440  float ipt[nmax - nmin];
441  float phiin[nmax - nmin];
442  float theta[nmax - nmin];
443 #pragma omp simd
444  for (int n = nmin; n < nmax; ++n) {
445  // if (std::abs(r-r0)<0.0001f) {
446  // dprint("distance less than 1mum, skip");
447  // continue;
448  // }
449 
450  xin[n - nmin] = inPar(n, 0, 0);
451  yin[n - nmin] = inPar(n, 1, 0);
452  ipt[n - nmin] = inPar(n, 3, 0);
453  phiin[n - nmin] = inPar(n, 4, 0);
454  theta[n - nmin] = inPar(n, 5, 0);
455 
456  //dprint(std::endl);
457  }
458 
459  //debug = true;
460  for (int n = nmin; n < nmax; ++n) {
461  dprint_np(n,
462  "input parameters"
463  << " inPar(n, 0, 0)=" << std::setprecision(9) << inPar(n, 0, 0) << " inPar(n, 1, 0)="
464  << std::setprecision(9) << inPar(n, 1, 0) << " inPar(n, 2, 0)=" << std::setprecision(9)
465  << inPar(n, 2, 0) << " inPar(n, 3, 0)=" << std::setprecision(9) << inPar(n, 3, 0)
466  << " inPar(n, 4, 0)=" << std::setprecision(9) << inPar(n, 4, 0)
467  << " inPar(n, 5, 0)=" << std::setprecision(9) << inPar(n, 5, 0));
468  }
469 
470  float kinv[nmax - nmin];
471  float pt[nmax - nmin];
472 #pragma omp simd
473  for (int n = nmin; n < nmax; ++n) {
474  kinv[n - nmin] = 1.f / k[n - nmin];
475  pt[n - nmin] = 1.f / ipt[n - nmin];
476  }
477  float D[nmax - nmin];
478  float cosa[nmax - nmin];
479  float sina[nmax - nmin];
480  float cosah[nmax - nmin];
481  float sinah[nmax - nmin];
482  float id[nmax - nmin];
483 
484 #pragma omp simd
485  for (int n = nmin; n < nmax; ++n) {
486  D[n - nmin] = 0.;
487  }
488 
489  //no trig approx here, phi can be large
490  float cosPorT[nmax - nmin];
491  float sinPorT[nmax - nmin];
492 #pragma omp simd
493  for (int n = nmin; n < nmax; ++n) {
494  cosPorT[n - nmin] = std::cos(phiin[n - nmin]);
495  }
496 #pragma omp simd
497  for (int n = nmin; n < nmax; ++n) {
498  sinPorT[n - nmin] = std::sin(phiin[n - nmin]);
499  }
500 
501  float pxin[nmax - nmin];
502  float pyin[nmax - nmin];
503 #pragma omp simd
504  for (int n = nmin; n < nmax; ++n) {
505  pxin[n - nmin] = cosPorT[n - nmin] * pt[n - nmin];
506  pyin[n - nmin] = sinPorT[n - nmin] * pt[n - nmin];
507  }
508 
509  for (int n = nmin; n < nmax; ++n) {
510  dprint_np(n,
511  "k=" << std::setprecision(9) << k[n - nmin] << " pxin=" << std::setprecision(9) << pxin[n - nmin]
512  << " pyin=" << std::setprecision(9) << pyin[n - nmin] << " cosPorT=" << std::setprecision(9)
513  << cosPorT[n - nmin] << " sinPorT=" << std::setprecision(9) << sinPorT[n - nmin]
514  << " pt=" << std::setprecision(9) << pt[n - nmin]);
515  }
516 
517  float dDdx[nmax - nmin];
518  float dDdy[nmax - nmin];
519  float dDdipt[nmax - nmin];
520  float dDdphi[nmax - nmin];
521 
522 #pragma omp simd
523  for (int n = nmin; n < nmax; ++n) {
524  dDdipt[n - nmin] = 0.;
525  dDdphi[n - nmin] = 0.;
526  }
527 #pragma omp simd
528  for (int n = nmin; n < nmax; ++n) {
529  //derivatives initialized to value for first iteration, i.e. distance = r-r0in
530  dDdx[n - nmin] = r0[n - nmin] > 0.f ? -xin[n - nmin] / r0[n - nmin] : 0.f;
531  }
532 
533 #pragma omp simd
534  for (int n = nmin; n < nmax; ++n) {
535  dDdy[n - nmin] = r0[n - nmin] > 0.f ? -yin[n - nmin] / r0[n - nmin] : 0.f;
536  }
537 
538  float oodotp[nmax - nmin];
539  float x[nmax - nmin];
540  float y[nmax - nmin];
541  float oor0[nmax - nmin];
542  float dadipt[nmax - nmin];
543  float dadx[nmax - nmin];
544  float dady[nmax - nmin];
545  float pxca[nmax - nmin];
546  float pxsa[nmax - nmin];
547  float pyca[nmax - nmin];
548  float pysa[nmax - nmin];
549  float tmp[nmax - nmin];
550  float tmpx[nmax - nmin];
551  float tmpy[nmax - nmin];
552  float pxinold[nmax - nmin];
553 
554  CMS_UNROLL_LOOP_COUNT(Config::Niter)
555  for (int i = 0; i < Config::Niter; ++i) {
556 #pragma omp simd
557  for (int n = nmin; n < nmax; ++n) {
558  //compute distance and path for the current iteration
559  r0[n - nmin] = hipo(outPar(n, 0, 0), outPar(n, 1, 0));
560  }
561 
562  // Use one over dot product of transverse momentum and radial
563  // direction to scale the step. Propagation is prevented from reaching
564  // too close to the apex (dotp > 0.2).
565  // - Can / should we come up with a better approximation?
566  // - Can / should take +/- curvature into account?
567 
568 #pragma omp simd
569  for (int n = nmin; n < nmax; ++n) {
570  oodotp[n - nmin] =
571  r0[n - nmin] * pt[n - nmin] / (pxin[n - nmin] * outPar(n, 0, 0) + pyin[n - nmin] * outPar(n, 1, 0));
572  }
573 
574 #pragma omp simd
575  for (int n = nmin; n < nmax; ++n) {
576  if (oodotp[n - nmin] > 5.0f || oodotp[n - nmin] < 0) // 0.2 is 78.5 deg
577  {
578  outFailFlag(n, 0, 0) = 1;
579  oodotp[n - nmin] = 0.0f;
580  } else if (r[n - nmin] - r0[n - nmin] < 0.0f && pt[n - nmin] < 1.0f) {
581  // Scale down the correction for low-pT ingoing tracks.
582  oodotp[n - nmin] = 1.0f + (oodotp[n - nmin] - 1.0f) * pt[n - nmin];
583  }
584  }
585 
586 #pragma omp simd
587  for (int n = nmin; n < nmax; ++n) {
588  // Can we come up with a better approximation?
589  // Should take +/- curvature into account.
590  id[n - nmin] = (r[n - nmin] - r0[n - nmin]) * oodotp[n - nmin];
591  }
592 
593 #pragma omp simd
594  for (int n = nmin; n < nmax; ++n) {
595  D[n - nmin] += id[n - nmin];
596  }
597 
598  if constexpr (Config::useTrigApprox) {
599 #if !defined(__INTEL_COMPILER)
600 #pragma omp simd
601 #endif
602  for (int n = nmin; n < nmax; ++n) {
603  sincos4(id[n - nmin] * ipt[n - nmin] * kinv[n - nmin] * 0.5f, sinah[n - nmin], cosah[n - nmin]);
604  }
605  } else {
606 #if !defined(__INTEL_COMPILER)
607 #pragma omp simd
608 #endif
609  for (int n = nmin; n < nmax; ++n) {
610  cosah[n - nmin] = std::cos(id[n - nmin] * ipt[n - nmin] * kinv[n - nmin] * 0.5f);
611  sinah[n - nmin] = std::sin(id[n - nmin] * ipt[n - nmin] * kinv[n - nmin] * 0.5f);
612  }
613  }
614 
615 #pragma omp simd
616  for (int n = nmin; n < nmax; ++n) {
617  cosa[n - nmin] = 1.f - 2.f * sinah[n - nmin] * sinah[n - nmin];
618  sina[n - nmin] = 2.f * sinah[n - nmin] * cosah[n - nmin];
619  }
620 
621  for (int n = nmin; n < nmax; ++n) {
622  dprint_np(n,
623  "Attempt propagation from r="
624  << r0[n - nmin] << " to r=" << r[n - nmin] << std::endl
625  << " x=" << xin[n - nmin] << " y=" << yin[n - nmin] << " z=" << inPar(n, 2, 0)
626  << " px=" << pxin[n - nmin] << " py=" << pyin[n - nmin]
627  << " pz=" << pt[n - nmin] * std::tan(theta[n - nmin]) << " q=" << inChg(n, 0, 0) << std::endl
628  << " r=" << std::setprecision(9) << r[n - nmin] << " r0=" << std::setprecision(9)
629  << r0[n - nmin] << " id=" << std::setprecision(9) << id[n - nmin]
630  << " dr=" << std::setprecision(9) << r[n - nmin] - r0[n - nmin] << " cosa=" << cosa[n - nmin]
631  << " sina=" << sina[n - nmin] << " dir_cos(rad,pT)=" << 1.0f / oodotp[n - nmin]);
632  }
633 
634  //update derivatives on total distance
635  if (i + 1 != Config::Niter) {
636 #pragma omp simd
637  for (int n = nmin; n < nmax; ++n) {
638  x[n - nmin] = outPar(n, 0, 0);
639  y[n - nmin] = outPar(n, 1, 0);
640  }
641 #pragma omp simd
642  for (int n = nmin; n < nmax; ++n) {
643  oor0[n - nmin] =
644  (r0[n - nmin] > 0.f && std::abs(r[n - nmin] - r0[n - nmin]) > 0.0001f) ? 1.f / r0[n - nmin] : 0.f;
645  }
646 #pragma omp simd
647  for (int n = nmin; n < nmax; ++n) {
648  dadipt[n - nmin] = id[n - nmin] * kinv[n - nmin];
649  dadx[n - nmin] = -x[n - nmin] * ipt[n - nmin] * kinv[n - nmin] * oor0[n - nmin];
650  dady[n - nmin] = -y[n - nmin] * ipt[n - nmin] * kinv[n - nmin] * oor0[n - nmin];
651  pxca[n - nmin] = pxin[n - nmin] * cosa[n - nmin];
652  pxsa[n - nmin] = pxin[n - nmin] * sina[n - nmin];
653  pyca[n - nmin] = pyin[n - nmin] * cosa[n - nmin];
654  pysa[n - nmin] = pyin[n - nmin] * sina[n - nmin];
655  tmpx[n - nmin] = k[n - nmin] * dadx[n - nmin];
656  }
657 
658 #pragma omp simd
659  for (int n = nmin; n < nmax; ++n) {
660  dDdx[n - nmin] -= (x[n - nmin] * (1.f + tmpx[n - nmin] * (pxca[n - nmin] - pysa[n - nmin])) +
661  y[n - nmin] * tmpx[n - nmin] * (pyca[n - nmin] + pxsa[n - nmin])) *
662  oor0[n - nmin];
663  }
664 
665 #pragma omp simd
666  for (int n = nmin; n < nmax; ++n) {
667  tmpy[n - nmin] = k[n - nmin] * dady[n - nmin];
668  }
669 #pragma omp simd
670  for (int n = nmin; n < nmax; ++n) {
671  dDdy[n - nmin] -= (x[n - nmin] * tmpy[n - nmin] * (pxca[n - nmin] - pysa[n - nmin]) +
672  y[n - nmin] * (1.f + tmpy[n - nmin] * (pyca[n - nmin] + pxsa[n - nmin]))) *
673  oor0[n - nmin];
674  }
675 #pragma omp simd
676  for (int n = nmin; n < nmax; ++n) {
677  //now r0 depends on ipt and phi as well
678  tmp[n - nmin] = dadipt[n - nmin] * ipt[n - nmin];
679  }
680 #pragma omp simd
681  for (int n = nmin; n < nmax; ++n) {
682  dDdipt[n - nmin] -= k[n - nmin] *
683  (x[n - nmin] * (pxca[n - nmin] * tmp[n - nmin] - pysa[n - nmin] * tmp[n - nmin] -
684  pyca[n - nmin] - pxsa[n - nmin] + pyin[n - nmin]) +
685  y[n - nmin] * (pyca[n - nmin] * tmp[n - nmin] + pxsa[n - nmin] * tmp[n - nmin] -
686  pysa[n - nmin] + pxca[n - nmin] - pxin[n - nmin])) *
687  pt[n - nmin] * oor0[n - nmin];
688  }
689 #pragma omp simd
690  for (int n = nmin; n < nmax; ++n) {
691  dDdphi[n - nmin] += k[n - nmin] *
692  (x[n - nmin] * (pysa[n - nmin] - pxin[n - nmin] + pxca[n - nmin]) -
693  y[n - nmin] * (pxsa[n - nmin] - pyin[n - nmin] + pyca[n - nmin])) *
694  oor0[n - nmin];
695  }
696  }
697 
698 #pragma omp simd
699  for (int n = nmin; n < nmax; ++n) {
700  //update parameters
701  outPar(n, 0, 0) = outPar(n, 0, 0) + 2.f * k[n - nmin] * sinah[n - nmin] *
702  (pxin[n - nmin] * cosah[n - nmin] - pyin[n - nmin] * sinah[n - nmin]);
703  outPar(n, 1, 0) = outPar(n, 1, 0) + 2.f * k[n - nmin] * sinah[n - nmin] *
704  (pyin[n - nmin] * cosah[n - nmin] + pxin[n - nmin] * sinah[n - nmin]);
705  pxinold[n - nmin] = pxin[n - nmin]; //copy before overwriting
706  pxin[n - nmin] = pxin[n - nmin] * cosa[n - nmin] - pyin[n - nmin] * sina[n - nmin];
707  pyin[n - nmin] = pyin[n - nmin] * cosa[n - nmin] + pxinold[n - nmin] * sina[n - nmin];
708  }
709  for (int n = nmin; n < nmax; ++n) {
710  dprint_np(n,
711  "outPar(n, 0, 0)=" << outPar(n, 0, 0) << " outPar(n, 1, 0)=" << outPar(n, 1, 0)
712  << " pxin=" << pxin[n - nmin] << " pyin=" << pyin[n - nmin]);
713  }
714  } // iteration loop
715 
716  float alpha[nmax - nmin];
717  float dadphi[nmax - nmin];
718 
719 #pragma omp simd
720  for (int n = nmin; n < nmax; ++n) {
721  alpha[n - nmin] = D[n - nmin] * ipt[n - nmin] * kinv[n - nmin];
722  dadx[n - nmin] = dDdx[n - nmin] * ipt[n - nmin] * kinv[n - nmin];
723  dady[n - nmin] = dDdy[n - nmin] * ipt[n - nmin] * kinv[n - nmin];
724  dadipt[n - nmin] = (ipt[n - nmin] * dDdipt[n - nmin] + D[n - nmin]) * kinv[n - nmin];
725  dadphi[n - nmin] = dDdphi[n - nmin] * ipt[n - nmin] * kinv[n - nmin];
726  }
727 
728  if constexpr (Config::useTrigApprox) {
729 #pragma omp simd
730  for (int n = nmin; n < nmax; ++n) {
731  sincos4(alpha[n - nmin], sina[n - nmin], cosa[n - nmin]);
732  }
733  } else {
734 #pragma omp simd
735  for (int n = nmin; n < nmax; ++n) {
736  cosa[n - nmin] = std::cos(alpha[n - nmin]);
737  sina[n - nmin] = std::sin(alpha[n - nmin]);
738  }
739  }
740 #pragma omp simd
741  for (int n = nmin; n < nmax; ++n) {
742  errorProp(n, 0, 0) = 1.f + k[n - nmin] * dadx[n - nmin] *
743  (cosPorT[n - nmin] * cosa[n - nmin] - sinPorT[n - nmin] * sina[n - nmin]) *
744  pt[n - nmin];
745  errorProp(n, 0, 1) = k[n - nmin] * dady[n - nmin] *
746  (cosPorT[n - nmin] * cosa[n - nmin] - sinPorT[n - nmin] * sina[n - nmin]) * pt[n - nmin];
747  errorProp(n, 0, 2) = 0.f;
748  errorProp(n, 0, 3) =
749  k[n - nmin] *
750  (cosPorT[n - nmin] * (ipt[n - nmin] * dadipt[n - nmin] * cosa[n - nmin] - sina[n - nmin]) +
751  sinPorT[n - nmin] * ((1.f - cosa[n - nmin]) - ipt[n - nmin] * dadipt[n - nmin] * sina[n - nmin])) *
752  pt[n - nmin] * pt[n - nmin];
753  errorProp(n, 0, 4) = k[n - nmin] *
754  (cosPorT[n - nmin] * dadphi[n - nmin] * cosa[n - nmin] -
755  sinPorT[n - nmin] * dadphi[n - nmin] * sina[n - nmin] - sinPorT[n - nmin] * sina[n - nmin] +
756  cosPorT[n - nmin] * cosa[n - nmin] - cosPorT[n - nmin]) *
757  pt[n - nmin];
758  errorProp(n, 0, 5) = 0.f;
759  }
760 
761 #pragma omp simd
762  for (int n = nmin; n < nmax; ++n) {
763  errorProp(n, 1, 0) = k[n - nmin] * dadx[n - nmin] *
764  (sinPorT[n - nmin] * cosa[n - nmin] + cosPorT[n - nmin] * sina[n - nmin]) * pt[n - nmin];
765  errorProp(n, 1, 1) = 1.f + k[n - nmin] * dady[n - nmin] *
766  (sinPorT[n - nmin] * cosa[n - nmin] + cosPorT[n - nmin] * sina[n - nmin]) *
767  pt[n - nmin];
768  errorProp(n, 1, 2) = 0.f;
769  errorProp(n, 1, 3) =
770  k[n - nmin] *
771  (sinPorT[n - nmin] * (ipt[n - nmin] * dadipt[n - nmin] * cosa[n - nmin] - sina[n - nmin]) +
772  cosPorT[n - nmin] * (ipt[n - nmin] * dadipt[n - nmin] * sina[n - nmin] - (1.f - cosa[n - nmin]))) *
773  pt[n - nmin] * pt[n - nmin];
774  errorProp(n, 1, 4) = k[n - nmin] *
775  (sinPorT[n - nmin] * dadphi[n - nmin] * cosa[n - nmin] +
776  cosPorT[n - nmin] * dadphi[n - nmin] * sina[n - nmin] + sinPorT[n - nmin] * cosa[n - nmin] +
777  cosPorT[n - nmin] * sina[n - nmin] - sinPorT[n - nmin]) *
778  pt[n - nmin];
779  errorProp(n, 1, 5) = 0.f;
780  }
781 
782 #pragma omp simd
783  for (int n = nmin; n < nmax; ++n) {
784  //no trig approx here, theta can be large
785  cosPorT[n - nmin] = std::cos(theta[n - nmin]);
786  }
787 
788 #pragma omp simd
789  for (int n = nmin; n < nmax; ++n) {
790  sinPorT[n - nmin] = std::sin(theta[n - nmin]);
791  }
792 
793 #pragma omp simd
794  for (int n = nmin; n < nmax; ++n) {
795  //redefine sinPorT as 1./sinPorT to reduce the number of temporaries
796  sinPorT[n - nmin] = 1.f / sinPorT[n - nmin];
797  }
798 
799 #pragma omp simd
800  for (int n = nmin; n < nmax; ++n) {
801  outPar(n, 2, 0) =
802  inPar(n, 2, 0) + k[n - nmin] * alpha[n - nmin] * cosPorT[n - nmin] * pt[n - nmin] * sinPorT[n - nmin];
803  errorProp(n, 2, 0) = k[n - nmin] * cosPorT[n - nmin] * dadx[n - nmin] * pt[n - nmin] * sinPorT[n - nmin];
804  errorProp(n, 2, 1) = k[n - nmin] * cosPorT[n - nmin] * dady[n - nmin] * pt[n - nmin] * sinPorT[n - nmin];
805  errorProp(n, 2, 2) = 1.f;
806  errorProp(n, 2, 3) = k[n - nmin] * cosPorT[n - nmin] * (ipt[n - nmin] * dadipt[n - nmin] - alpha[n - nmin]) *
807  pt[n - nmin] * pt[n - nmin] * sinPorT[n - nmin];
808  errorProp(n, 2, 4) = k[n - nmin] * dadphi[n - nmin] * cosPorT[n - nmin] * pt[n - nmin] * sinPorT[n - nmin];
809  errorProp(n, 2, 5) = -k[n - nmin] * alpha[n - nmin] * pt[n - nmin] * sinPorT[n - nmin] * sinPorT[n - nmin];
810  }
811 
812 #pragma omp simd
813  for (int n = nmin; n < nmax; ++n) {
814  outPar(n, 3, 0) = ipt[n - nmin];
815  errorProp(n, 3, 0) = 0.f;
816  errorProp(n, 3, 1) = 0.f;
817  errorProp(n, 3, 2) = 0.f;
818  errorProp(n, 3, 3) = 1.f;
819  errorProp(n, 3, 4) = 0.f;
820  errorProp(n, 3, 5) = 0.f;
821  }
822 
823 #pragma omp simd
824  for (int n = nmin; n < nmax; ++n) {
825  outPar(n, 4, 0) = inPar(n, 4, 0) + alpha[n - nmin];
826  errorProp(n, 4, 0) = dadx[n - nmin];
827  errorProp(n, 4, 1) = dady[n - nmin];
828  errorProp(n, 4, 2) = 0.f;
829  errorProp(n, 4, 3) = dadipt[n - nmin];
830  errorProp(n, 4, 4) = 1.f + dadphi[n - nmin];
831  errorProp(n, 4, 5) = 0.f;
832  }
833 
834 #pragma omp simd
835  for (int n = nmin; n < nmax; ++n) {
836  outPar(n, 5, 0) = theta[n - nmin];
837  errorProp(n, 5, 0) = 0.f;
838  errorProp(n, 5, 1) = 0.f;
839  errorProp(n, 5, 2) = 0.f;
840  errorProp(n, 5, 3) = 0.f;
841  errorProp(n, 5, 4) = 0.f;
842  errorProp(n, 5, 5) = 1.f;
843  }
844 
845  for (int n = nmin; n < nmax; ++n) {
846  dprint_np(
847  n,
848  "propagation end, dump parameters\n"
849  << " D = " << D[n - nmin] << " alpha = " << alpha[n - nmin] << " kinv = " << kinv[n - nmin] << std::endl
850  << " pos = " << outPar(n, 0, 0) << " " << outPar(n, 1, 0) << " " << outPar(n, 2, 0) << "\t\t r="
851  << std::sqrt(outPar(n, 0, 0) * outPar(n, 0, 0) + outPar(n, 1, 0) * outPar(n, 1, 0)) << std::endl
852  << " mom = " << outPar(n, 3, 0) << " " << outPar(n, 4, 0) << " " << outPar(n, 5, 0) << std::endl
853  << " cart= " << std::cos(outPar(n, 4, 0)) / outPar(n, 3, 0) << " "
854  << std::sin(outPar(n, 4, 0)) / outPar(n, 3, 0) << " " << 1. / (outPar(n, 3, 0) * tan(outPar(n, 5, 0)))
855  << "\t\tpT=" << 1. / std::abs(outPar(n, 3, 0)) << std::endl);
856  }
857 
858 #ifdef DEBUG
859  for (int n = nmin; n < nmax; ++n) {
860  if (debug && g_debug && n < N_proc) {
861  dmutex_guard;
862  std::cout << n << ": jacobian" << std::endl;
863  printf("%5f %5f %5f %5f %5f %5f\n",
864  errorProp(n, 0, 0),
865  errorProp(n, 0, 1),
866  errorProp(n, 0, 2),
867  errorProp(n, 0, 3),
868  errorProp(n, 0, 4),
869  errorProp(n, 0, 5));
870  printf("%5f %5f %5f %5f %5f %5f\n",
871  errorProp(n, 1, 0),
872  errorProp(n, 1, 1),
873  errorProp(n, 1, 2),
874  errorProp(n, 1, 3),
875  errorProp(n, 1, 4),
876  errorProp(n, 1, 5));
877  printf("%5f %5f %5f %5f %5f %5f\n",
878  errorProp(n, 2, 0),
879  errorProp(n, 2, 1),
880  errorProp(n, 2, 2),
881  errorProp(n, 2, 3),
882  errorProp(n, 2, 4),
883  errorProp(n, 2, 5));
884  printf("%5f %5f %5f %5f %5f %5f\n",
885  errorProp(n, 3, 0),
886  errorProp(n, 3, 1),
887  errorProp(n, 3, 2),
888  errorProp(n, 3, 3),
889  errorProp(n, 3, 4),
890  errorProp(n, 3, 5));
891  printf("%5f %5f %5f %5f %5f %5f\n",
892  errorProp(n, 4, 0),
893  errorProp(n, 4, 1),
894  errorProp(n, 4, 2),
895  errorProp(n, 4, 3),
896  errorProp(n, 4, 4),
897  errorProp(n, 4, 5));
898  printf("%5f %5f %5f %5f %5f %5f\n",
899  errorProp(n, 5, 0),
900  errorProp(n, 5, 1),
901  errorProp(n, 5, 2),
902  errorProp(n, 5, 3),
903  errorProp(n, 5, 4),
904  errorProp(n, 5, 5));
905  printf("\n");
906  }
907  }
908 #endif
909  }
910 
911 } // namespace
912 
913 // END STUFF FROM PropagationMPlex.icc
914 // ============================================================================
915 
916 namespace mkfit {
917 
918  void helixAtRFromIterativeCCS(const MPlexLV& inPar,
919  const MPlexQI& inChg,
920  const MPlexQF& msRad,
921  MPlexLV& outPar,
922  MPlexLL& errorProp,
923  MPlexQI& outFailFlag,
924  const int N_proc,
925  const PropagationFlags& pflags) {
926  errorProp.setVal(0.f);
927  outFailFlag.setVal(0.f);
928 
929  helixAtRFromIterativeCCS_impl(inPar, inChg, msRad, outPar, errorProp, outFailFlag, 0, NN, N_proc, pflags);
930  }
931 
932  void propagateHelixToRMPlex(const MPlexLS& inErr,
933  const MPlexLV& inPar,
934  const MPlexQI& inChg,
935  const MPlexQF& msRad,
936  MPlexLS& outErr,
937  MPlexLV& outPar,
938  MPlexQI& outFailFlag,
939  const int N_proc,
940  const PropagationFlags& pflags,
941  const MPlexQI* noMatEffPtr) {
942  // bool debug = true;
943 
944  // This is used further down when calculating similarity with errorProp (and before in DEBUG).
945  // MT: I don't think this really needed if we use inErr where required.
946  outErr = inErr;
947  // This requirement for helixAtRFromIterativeCCS_impl() and for helixAtRFromIterativeCCSFullJac().
948  // MT: This should be properly handled in both functions (expecting input in output parameters sucks).
949  outPar = inPar;
950 
951  MPlexLL errorProp;
952 
953  helixAtRFromIterativeCCS(inPar, inChg, msRad, outPar, errorProp, outFailFlag, N_proc, pflags);
954 
955 #ifdef DEBUG
956  if (debug && g_debug) {
957  for (int kk = 0; kk < N_proc; ++kk) {
958  dprintf("outErr before prop %d\n", kk);
959  for (int i = 0; i < 6; ++i) {
960  for (int j = 0; j < 6; ++j)
961  dprintf("%8f ", outErr.At(kk, i, j));
962  dprintf("\n");
963  }
964  dprintf("\n");
965 
966  dprintf("errorProp %d\n", kk);
967  for (int i = 0; i < 6; ++i) {
968  for (int j = 0; j < 6; ++j)
969  dprintf("%8f ", errorProp.At(kk, i, j));
970  dprintf("\n");
971  }
972  dprintf("\n");
973  }
974  }
975 #endif
976 
977  // MultHelixProp can be optimized for CCS coordinates, see GenMPlexOps.pl
978  MPlexLL temp;
979  MultHelixProp(errorProp, outErr, temp);
980  MultHelixPropTransp(errorProp, temp, outErr);
981  // can replace with: MultHelixPropFull(errorProp, outErr, temp); MultHelixPropTranspFull(errorProp, temp, outErr);
982 
983 #ifdef DEBUG
984  if (debug && g_debug) {
985  for (int kk = 0; kk < N_proc; ++kk) {
986  dprintf("outErr %d\n", kk);
987  for (int i = 0; i < 6; ++i) {
988  for (int j = 0; j < 6; ++j)
989  dprintf("%8f ", outErr.constAt(kk, i, j));
990  dprintf("\n");
991  }
992  dprintf("\n");
993  }
994  }
995 #endif
996 
997  if (pflags.apply_material) {
998  MPlexQF hitsRl;
999  MPlexQF hitsXi;
1000  MPlexQF propSign;
1001 
1002  const TrackerInfo& tinfo = *pflags.tracker_info;
1003 
1004 #pragma omp simd
1005  for (int n = 0; n < NN; ++n) {
1006  if (n < N_proc) {
1007  if (outFailFlag(n, 0, 0) || (noMatEffPtr && noMatEffPtr->constAt(n, 0, 0))) {
1008  hitsRl(n, 0, 0) = 0.f;
1009  hitsXi(n, 0, 0) = 0.f;
1010  } else {
1011  auto mat = tinfo.material_checked(std::abs(outPar(n, 2, 0)), msRad(n, 0, 0));
1012  hitsRl(n, 0, 0) = mat.radl;
1013  hitsXi(n, 0, 0) = mat.bbxi;
1014  }
1015  const float r0 = hipo(inPar(n, 0, 0), inPar(n, 1, 0));
1016  const float r = msRad(n, 0, 0);
1017  propSign(n, 0, 0) = (r > r0 ? 1. : -1.);
1018  }
1019  }
1020  MPlexHV plNrm;
1021 #pragma omp simd
1022  for (int n = 0; n < NN; ++n) {
1023  plNrm(n, 0, 0) = std::cos(outPar.constAt(n, 4, 0));
1024  plNrm(n, 1, 0) = std::sin(outPar.constAt(n, 4, 0));
1025  plNrm(n, 2, 0) = 0.f;
1026  }
1027  applyMaterialEffects(hitsRl, hitsXi, propSign, plNrm, outErr, outPar, N_proc);
1028 #ifdef DEBUG
1029  if (debug && g_debug) {
1030  for (int kk = 0; kk < N_proc; ++kk) {
1031  dprintf("propSign %d\n", kk);
1032  for (int i = 0; i < 1; ++i) {
1033  dprintf("%8f ", propSign.constAt(kk, i, 0));
1034  }
1035  dprintf("\n");
1036  dprintf("plNrm %d\n", kk);
1037  for (int i = 0; i < 3; ++i) {
1038  dprintf("%8f ", plNrm.constAt(kk, i, 0));
1039  }
1040  dprintf("\n");
1041  dprintf("outErr(after material) %d\n", kk);
1042  for (int i = 0; i < 6; ++i) {
1043  for (int j = 0; j < 6; ++j)
1044  dprintf("%8f ", outErr.constAt(kk, i, j));
1045  dprintf("\n");
1046  }
1047  dprintf("\n");
1048  }
1049  }
1050 #endif
1051  }
1052 
1053  squashPhiMPlex(outPar, N_proc); // ensure phi is between |pi|
1054 
1055  // Matriplex version of:
1056  // result.errors = ROOT::Math::Similarity(errorProp, outErr);
1057 
1058  /*
1059  // To be used with: MPT_DIM = 1
1060  if (fabs(sqrt(outPar[0]*outPar[0]+outPar[1]*outPar[1]) - msRad[0]) > 0.0001)
1061  {
1062  std::cout << "DID NOT GET TO R, FailFlag=" << failFlag[0]
1063  << " dR=" << msRad[0] - std::hypot(outPar[0],outPar[1])
1064  << " r=" << msRad[0] << " rin=" << std::hypot(inPar[0],inPar[1]) << " rout=" << std::hypot(outPar[0],outPar[1])
1065  << std::endl;
1066  // std::cout << " pt=" << pt << " pz=" << inPar.At(n, 2) << std::endl;
1067  }
1068  */
1069 
1070  // PROP-FAIL-ENABLE To keep physics changes minimal, we always restore the
1071  // state to input when propagation fails -- as was the default before.
1072  // if (pflags.copy_input_state_on_fail) {
1073  for (int i = 0; i < N_proc; ++i) {
1074  if (outFailFlag(i, 0, 0)) {
1075  outPar.copySlot(i, inPar);
1076  outErr.copySlot(i, inErr);
1077  }
1078  }
1079  // }
1080  }
1081 
1082 } // end namespace mkfit
Matriplex::sincos4
void sincos4(const MPlex< T, D1, D2, N > &a, MPlex< T, D1, D2, N > &s, MPlex< T, D1, D2, N > &c)
Definition: Matriplex.h:696
mkfit::MPlexLL
Matriplex::Matriplex< float, LL, LL, NN > MPlexLL
Definition: Matrix.h:53
mps_fire.i
i
Definition: mps_fire.py:429
mkfit::propagateHelixToRMPlex
void propagateHelixToRMPlex(const MPlexLS &inErr, const MPlexLV &inPar, const MPlexQI &inChg, const MPlexQF &msRad, MPlexLS &outErr, MPlexLV &outPar, MPlexQI &outFailFlag, const int N_proc, const PropagationFlags &pflags, const MPlexQI *noMatEffPtr)
Definition: PropagationMPlex.cc:932
B
Definition: APVGainStruct.h:7
dprint_np
#define dprint_np(n, x)
Definition: Debug.h:96
tauSpinnerTable_cfi.theta
theta
Definition: tauSpinnerTable_cfi.py:7
mkfit::MPlexHV
Matriplex::Matriplex< float, HH, 1, NN > MPlexHV
Definition: Matrix.h:58
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
CMS_UNROLL_LOOP_COUNT
#define CMS_UNROLL_LOOP_COUNT(N)
Definition: CMSUnrollLoop.h:48
mkfit::propagateLineToRMPlex
void propagateLineToRMPlex(const MPlexLS &psErr, const MPlexLV &psPar, const MPlexHS &msErr, const MPlexHV &msPar, MPlexLS &outErr, MPlexLV &outPar, const int N_proc)
Definition: PropagationMPlex.cc:17
mkfit::PropagationFlags::tracker_info
const TrackerInfo * tracker_info
Definition: PropagationConfig.h:17
mkfit::applyMaterialEffects
void applyMaterialEffects(const MPlexQF &hitsRl, const MPlexQF &hitsXi, const MPlexQF &propSign, const MPlexHV &plNrm, MPlexLS &outErr, MPlexLV &outPar, const int N_proc)
Definition: PropagationMPlexCommon.cc:77
ALPAKA_ACCELERATOR_NAMESPACE::brokenline::constexpr
if constexpr(n > 3)
Definition: BrokenLine.h:164
mkfit::MPlexLV
Matriplex::Matriplex< float, LL, 1, NN > MPlexLV
Definition: Matrix.h:54
mkfit::squashPhiMPlex
void squashPhiMPlex(MPlexLV &par, const int N_proc)
Definition: PropagationMPlex.h:10
DiDispStaMuonMonitor_cfi.pt
pt
Definition: DiDispStaMuonMonitor_cfi.py:39
packedPFCandidateRefMixer_cfi.pf
pf
Definition: packedPFCandidateRefMixer_cfi.py:4
mkfit::g_debug
bool g_debug
Definition: Debug.cc:2
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:23
mkfit::NN
constexpr Matriplex::idx_t NN
Definition: Matrix.h:48
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
mkfit::Config::Bfield
constexpr float Bfield
Definition: Config.h:60
mkfit::helixAtRFromIterativeCCSFullJac
void helixAtRFromIterativeCCSFullJac(const MPlexLV &inPar, const MPlexQI &inChg, const MPlexQF &msRad, MPlexLV &outPar, MPlexLL &errorProp, const int N_proc)
Definition: PropagationMPlex.cc:201
funct::tan
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
f
double f[11][100]
Definition: MuScleFitUtils.cc:78
mkfit::Config::useTrigApprox
constexpr bool useTrigApprox
Definition: Config.h:51
relval_data_highstats.y
y
Definition: relval_data_highstats.py:28
correctionTermsCaloMet_cff.C
C
Definition: correctionTermsCaloMet_cff.py:34
mkfit::MPlexQI
Matriplex::Matriplex< int, 1, 1, NN > MPlexQI
Definition: Matrix.h:81
debug
#define debug
Definition: HDRShower.cc:19
mkfit::Const::sol
constexpr float sol
Definition: Config.h:13
mkfit::MPlexQF
Matriplex::Matriplex< float, 1, 1, NN > MPlexQF
Definition: Matrix.h:80
N
#define N
Definition: blowfish.cc:9
funct::D
DecomposeProduct< arg, typename Div::arg > D
Definition: Factorize.h:141
mkfit::hipo
float hipo(float x, float y)
Definition: Matrix.h:9
mkfit::Config::bFieldFromZR
float bFieldFromZR(const float z, const float r)
Definition: Config.h:131
create_idmaps.n
n
Definition: create_idmaps.py:21
b
double b
Definition: hdecay.h:120
mkfit::MPlexHS
Matriplex::MatriplexSym< float, HH, NN > MPlexHS
Definition: Matrix.h:59
GetRecoTauVFromDQM_MC_cff.kk
kk
Definition: GetRecoTauVFromDQM_MC_cff.py:84
TrackingDataMCValidation_Standalone_cff.nmin
nmin
Definition: TrackingDataMCValidation_Standalone_cff.py:26
a
double a
Definition: hdecay.h:121
mkfit::MPlexLS
Matriplex::MatriplexSym< float, LL, NN > MPlexLS
Definition: Matrix.h:55
x
float x
Definition: beamSpotDipStandalone.cc:55
mkfit::TrackerInfo::material_checked
Material material_checked(float z, float r) const
Definition: TrackerInfo.h:275
gather_cfg.cout
cout
Definition: gather_cfg.py:144
A
Definition: APVGainStruct.h:7
createJobs.tmp
tmp
align.sh
Definition: createJobs.py:716
mkfit::helixAtRFromIterativeCCS
void helixAtRFromIterativeCCS(const MPlexLV &inPar, const MPlexQI &inChg, const MPlexQF &msRad, MPlexLV &outPar, MPlexLL &errorProp, MPlexQI &outFailFlag, const int N_proc, const PropagationFlags &pflags)
Definition: PropagationMPlex.cc:918
T
long double T
Definition: Basic3DVectorLD.h:48
mkfit::sincos4
void sincos4(const float x, float &sin, float &cos)
Definition: Matrix.h:13
ASSUME_ALIGNED
#define ASSUME_ALIGNED(a, b)
Definition: MatriplexCommon.h:94
dprintf
#define dprintf(...)
Definition: Debug.h:98
AlCaHLTBitMon_ParallelJobs.p
def p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
mkfit::Config::Niter
constexpr int Niter
Definition: Config.h:50
Generated for CMSSW Reference Manual by doxygen 1.8.14