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