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