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 (debug && g_debug && 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  MPlexQI& outFailFlag,
494  const int N_proc,
495  const PropagationFlags pflags,
496  const MPlexQI* noMatEffPtr) {
497  // bool debug = true;
498 
499  // This is used further down when calculating similarity with errorProp (and before in DEBUG).
500  // MT: I don't think this really needed if we use inErr where required.
501  outErr = inErr;
502  // This requirement for helixAtRFromIterativeCCS_impl() and for helixAtRFromIterativeCCSFullJac().
503  // MT: This should be properly handled in both functions (expecting input in output parameters sucks).
504  outPar = inPar;
505 
506  MPlexLL errorProp;
507 
508  helixAtRFromIterativeCCS(inPar, inChg, msRad, outPar, errorProp, outFailFlag, N_proc, pflags);
509 
510 #ifdef DEBUG
511  if (debug && g_debug) {
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  dprintf("\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  dprintf("\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 < NN; ++n) {
538  if (n >= N_proc || (outFailFlag(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  // PROP-FAIL-ENABLE To keep physics changes minimal, we always restore the
581  // state to input when propagation fails -- as was the default before.
582  // if (pflags.copy_input_state_on_fail) {
583  for (int i = 0; i < N_proc; ++i) {
584  if (outFailFlag(i, 0, 0)) {
585  outPar.copySlot(i, inPar);
586  outErr.copySlot(i, inErr);
587  }
588  }
589  // }
590  }
591 
592  //==============================================================================
593 
594  void propagateHelixToZMPlex(const MPlexLS& inErr,
595  const MPlexLV& inPar,
596  const MPlexQI& inChg,
597  const MPlexQF& msZ,
598  MPlexLS& outErr,
599  MPlexLV& outPar,
600  MPlexQI& outFailFlag,
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, outFailFlag, N_proc, pflags);
612 
613 #ifdef DEBUG
614  if (debug && g_debug) {
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  dprintf("\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  dprintf("\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 < NN; ++n) {
641  if (n >= N_proc || (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  // PROP-FAIL-ENABLE To keep physics changes minimal, we always restore the
670  // state to input when propagation fails -- as was the default before.
671  // if (pflags.copy_input_state_on_fail) {
672  for (int i = 0; i < N_proc; ++i) {
673  if (outFailFlag(i, 0, 0)) {
674  outPar.copySlot(i, inPar);
675  outErr.copySlot(i, inErr);
676  }
677  }
678  // }
679 
680  // This dump is now out of its place as similarity is done with matriplex ops.
681  /*
682 #ifdef DEBUG
683  {
684  dmutex_guard;
685  for (int kk = 0; kk < N_proc; ++kk)
686  {
687  dprintf("outErr %d\n", kk);
688  for (int i = 0; i < 6; ++i) { for (int j = 0; j < 6; ++j)
689  dprintf("%8f ", outErr.At(kk,i,j)); printf("\n");
690  } dprintf("\n");
691 
692  dprintf("outPar %d\n", kk);
693  for (int i = 0; i < 6; ++i) {
694  dprintf("%8f ", outPar.At(kk,i,0)); printf("\n");
695  } dprintf("\n");
696  if (std::abs(outPar.At(kk,2,0) - msZ.constAt(kk, 0, 0)) > 0.0001) {
697  float pt = 1.0f / inPar.constAt(kk,3,0);
698  dprint_np(kk, "DID NOT GET TO Z, dZ=" << std::abs(outPar.At(kk,2,0) - msZ.constAt(kk, 0, 0))
699  << " z=" << msZ.constAt(kk, 0, 0) << " zin=" << inPar.constAt(kk,2,0) << " zout=" << outPar.At(kk,2,0) << std::endl
700  << "pt=" << pt << " pz=" << pt/std::tan(inPar.constAt(kk,5,0)));
701  }
702  }
703  }
704 #endif
705  */
706  }
707 
708  void helixAtZ(const MPlexLV& inPar,
709  const MPlexQI& inChg,
710  const MPlexQF& msZ,
711  MPlexLV& outPar,
712  MPlexLL& errorProp,
713  MPlexQI& outFailFlag,
714  const int N_proc,
715  const PropagationFlags pflags) {
716  errorProp.setVal(0.f);
717 
718 #pragma omp simd
719  for (int n = 0; n < NN; ++n) {
720  //initialize erroProp to identity matrix, except element 2,2 which is zero
721  errorProp(n, 0, 0) = 1.f;
722  errorProp(n, 1, 1) = 1.f;
723  errorProp(n, 3, 3) = 1.f;
724  errorProp(n, 4, 4) = 1.f;
725  errorProp(n, 5, 5) = 1.f;
726  }
727  float zout[NN];
728  float zin[NN];
729  float ipt[NN];
730  float phiin[NN];
731  float theta[NN];
732 #pragma omp simd
733  for (int n = 0; n < NN; ++n) {
734  //initialize erroProp to identity matrix, except element 2,2 which is zero
735  zout[n] = msZ.constAt(n, 0, 0);
736  zin[n] = inPar.constAt(n, 2, 0);
737  ipt[n] = inPar.constAt(n, 3, 0);
738  phiin[n] = inPar.constAt(n, 4, 0);
739  theta[n] = inPar.constAt(n, 5, 0);
740  }
741 
742  float k[NN];
743  if (pflags.use_param_b_field) {
744 #pragma omp simd
745  for (int n = 0; n < NN; ++n) {
746  k[n] = inChg.constAt(n, 0, 0) * 100.f /
747  (-Const::sol * Config::bFieldFromZR(zin[n], hipo(inPar.constAt(n, 0, 0), inPar.constAt(n, 1, 0))));
748  }
749  } else {
750 #pragma omp simd
751  for (int n = 0; n < NN; ++n) {
752  k[n] = inChg.constAt(n, 0, 0) * 100.f / (-Const::sol * Config::Bfield);
753  }
754  }
755 
756  float kinv[NN];
757 #pragma omp simd
758  for (int n = 0; n < NN; ++n) {
759  kinv[n] = 1.f / k[n];
760  }
761 
762 #pragma omp simd
763  for (int n = 0; n < NN; ++n) {
764  dprint_np(n,
765  std::endl
766  << "input parameters"
767  << " inPar.constAt(n, 0, 0)=" << std::setprecision(9) << inPar.constAt(n, 0, 0)
768  << " inPar.constAt(n, 1, 0)=" << std::setprecision(9) << inPar.constAt(n, 1, 0)
769  << " inPar.constAt(n, 2, 0)=" << std::setprecision(9) << inPar.constAt(n, 2, 0)
770  << " inPar.constAt(n, 3, 0)=" << std::setprecision(9) << inPar.constAt(n, 3, 0)
771  << " inPar.constAt(n, 4, 0)=" << std::setprecision(9) << inPar.constAt(n, 4, 0)
772  << " inPar.constAt(n, 5, 0)=" << std::setprecision(9) << inPar.constAt(n, 5, 0));
773  }
774 
775  float pt[NN];
776 #pragma omp simd
777  for (int n = 0; n < NN; ++n) {
778  pt[n] = 1.f / ipt[n];
779  }
780 
781  //no trig approx here, phi can be large
782  float cosP[NN];
783  float sinP[NN];
784 #pragma omp simd
785  for (int n = 0; n < NN; ++n) {
786  cosP[n] = std::cos(phiin[n]);
787  }
788 
789 #pragma omp simd
790  for (int n = 0; n < NN; ++n) {
791  sinP[n] = std::sin(phiin[n]);
792  }
793 
794  float cosT[NN];
795  float sinT[NN];
796 #pragma omp simd
797  for (int n = 0; n < NN; ++n) {
798  cosT[n] = std::cos(theta[n]);
799  }
800 
801 #pragma omp simd
802  for (int n = 0; n < NN; ++n) {
803  sinT[n] = std::sin(theta[n]);
804  }
805 
806  float tanT[NN];
807  float icos2T[NN];
808  float pxin[NN];
809  float pyin[NN];
810 #pragma omp simd
811  for (int n = 0; n < NN; ++n) {
812  tanT[n] = sinT[n] / cosT[n];
813  icos2T[n] = 1.f / (cosT[n] * cosT[n]);
814  pxin[n] = cosP[n] * pt[n];
815  pyin[n] = sinP[n] * pt[n];
816  }
817 #pragma omp simd
818  for (int n = 0; n < NN; ++n) {
819  //fixme, make this printout useful for propagation to z
820  dprint_np(n,
821  std::endl
822  << "k=" << std::setprecision(9) << k[n] << " pxin=" << std::setprecision(9) << pxin[n]
823  << " pyin=" << std::setprecision(9) << pyin[n] << " cosP=" << std::setprecision(9) << cosP[n]
824  << " sinP=" << std::setprecision(9) << sinP[n] << " pt=" << std::setprecision(9) << pt[n]);
825  }
826  float deltaZ[NN];
827  float alpha[NN];
828 #pragma omp simd
829  for (int n = 0; n < NN; ++n) {
830  deltaZ[n] = zout[n] - zin[n];
831  alpha[n] = deltaZ[n] * tanT[n] * ipt[n] * kinv[n];
832  }
833 
834  float cosahTmp[NN];
835  float sinahTmp[NN];
836  if constexpr (Config::useTrigApprox) {
837 #if !defined(__INTEL_COMPILER)
838 #pragma omp simd
839 #endif
840  for (int n = 0; n < NN; ++n) {
841  sincos4(alpha[n] * 0.5f, sinahTmp[n], cosahTmp[n]);
842  }
843  } else {
844 #if !defined(__INTEL_COMPILER)
845 #pragma omp simd
846 #endif
847  for (int n = 0; n < NN; ++n) {
848  cosahTmp[n] = std::cos(alpha[n] * 0.5f);
849  }
850 #if !defined(__INTEL_COMPILER)
851 #pragma omp simd
852 #endif
853  for (int n = 0; n < NN; ++n) {
854  sinahTmp[n] = std::sin(alpha[n] * 0.5f);
855  }
856  }
857 
858  float cosah[NN];
859  float sinah[NN];
860  float cosa[NN];
861  float sina[NN];
862 #pragma omp simd
863  for (int n = 0; n < NN; ++n) {
864  cosah[n] = cosahTmp[n];
865  sinah[n] = sinahTmp[n];
866  cosa[n] = 1.f - 2.f * sinah[n] * sinah[n];
867  sina[n] = 2.f * sinah[n] * cosah[n];
868  }
869 //update parameters
870 #pragma omp simd
871  for (int n = 0; n < NN; ++n) {
872  outPar.At(n, 0, 0) = outPar.At(n, 0, 0) + 2.f * k[n] * sinah[n] * (pxin[n] * cosah[n] - pyin[n] * sinah[n]);
873  outPar.At(n, 1, 0) = outPar.At(n, 1, 0) + 2.f * k[n] * sinah[n] * (pyin[n] * cosah[n] + pxin[n] * sinah[n]);
874  outPar.At(n, 2, 0) = zout[n];
875  outPar.At(n, 4, 0) = phiin[n] + alpha[n];
876  }
877 
878 #pragma omp simd
879  for (int n = 0; n < NN; ++n) {
880  dprint_np(n,
881  std::endl
882  << "outPar.At(n, 0, 0)=" << outPar.At(n, 0, 0) << " outPar.At(n, 1, 0)=" << outPar.At(n, 1, 0)
883  << " pxin=" << pxin[n] << " pyin=" << pyin[n]);
884  }
885 
886  float pxcaMpysa[NN];
887 #pragma omp simd
888  for (int n = 0; n < NN; ++n) {
889  pxcaMpysa[n] = pxin[n] * cosa[n] - pyin[n] * sina[n];
890  }
891 
892 #pragma omp simd
893  for (int n = 0; n < NN; ++n) {
894  errorProp(n, 0, 2) = -tanT[n] * ipt[n] * pxcaMpysa[n];
895  errorProp(n, 0, 3) =
896  k[n] * pt[n] * pt[n] *
897  (cosP[n] * (alpha[n] * cosa[n] - sina[n]) + sinP[n] * 2.f * sinah[n] * (sinah[n] - alpha[n] * cosah[n]));
898  errorProp(n, 0, 4) = -2.f * k[n] * pt[n] * sinah[n] * (sinP[n] * cosah[n] + cosP[n] * sinah[n]);
899  errorProp(n, 0, 5) = deltaZ[n] * ipt[n] * pxcaMpysa[n] * icos2T[n];
900  }
901 
902  float pycaPpxsa[NN];
903 #pragma omp simd
904  for (int n = 0; n < NN; ++n) {
905  pycaPpxsa[n] = pyin[n] * cosa[n] + pxin[n] * sina[n];
906  }
907 
908 #pragma omp simd
909  for (int n = 0; n < NN; ++n) {
910  errorProp(n, 1, 2) = -tanT[n] * ipt[n] * pycaPpxsa[n];
911  errorProp(n, 1, 3) =
912  k[n] * pt[n] * pt[n] *
913  (sinP[n] * (alpha[n] * cosa[n] - sina[n]) - cosP[n] * 2.f * sinah[n] * (sinah[n] - alpha[n] * cosah[n]));
914  errorProp(n, 1, 4) = 2.f * k[n] * pt[n] * sinah[n] * (cosP[n] * cosah[n] - sinP[n] * sinah[n]);
915  errorProp(n, 1, 5) = deltaZ[n] * ipt[n] * pycaPpxsa[n] * icos2T[n];
916  }
917 
918 #pragma omp simd
919  for (int n = 0; n < NN; ++n) {
920  errorProp(n, 4, 2) = -ipt[n] * tanT[n] * kinv[n];
921  errorProp(n, 4, 3) = tanT[n] * deltaZ[n] * kinv[n];
922  errorProp(n, 4, 5) = ipt[n] * deltaZ[n] * kinv[n] * icos2T[n];
923  }
924 
925 #pragma omp simd
926  for (int n = 0; n < NN; ++n) {
927  dprint_np(
928  n,
929  "propagation end, dump parameters"
930  << std::endl
931  << "pos = " << outPar.At(n, 0, 0) << " " << outPar.At(n, 1, 0) << " " << outPar.At(n, 2, 0) << std::endl
932  << "mom = " << std::cos(outPar.At(n, 4, 0)) / outPar.At(n, 3, 0) << " "
933  << std::sin(outPar.At(n, 4, 0)) / outPar.At(n, 3, 0) << " "
934  << 1. / (outPar.At(n, 3, 0) * tan(outPar.At(n, 5, 0)))
935  << " r=" << std::sqrt(outPar.At(n, 0, 0) * outPar.At(n, 0, 0) + outPar.At(n, 1, 0) * outPar.At(n, 1, 0))
936  << " pT=" << 1. / std::abs(outPar.At(n, 3, 0)) << std::endl);
937  }
938 
939  // PROP-FAIL-ENABLE Disabled to keep physics changes minimal.
940  // To be reviewed, enabled and processed accordingly elsewhere.
941  /*
942  // Check for errors, set fail-flag.
943  for (int n = 0; n < NN; ++n) {
944  // We propagate for alpha: mark fail when prop angle more than pi/2
945  if (std::abs(alpha[n]) > 1.57) {
946  dprintf("helixAtZ: more than quarter turn, alpha = %f\n", alpha[n]);
947  outFailFlag[n] = 1;
948  } else {
949  // Have we reached desired z? We can't know, we copy desired z to actual z.
950  // Are we close to apex? Same condition as in propToR, 12.5 deg, cos(78.5deg) = 0.2
951  float dotp = (outPar.At(n, 0, 0) * std::cos(outPar.At(n, 4, 0)) +
952  outPar.At(n, 1, 0) * std::sin(outPar.At(n, 4, 0))) /
953  std::hypot(outPar.At(n, 0, 0), outPar.At(n, 1, 0));
954  if (dotp < 0.2 || dotp < 0) {
955  dprintf("helixAtZ: dot product bad, dotp = %f\n", dotp);
956  outFailFlag[n] = 1;
957  }
958  }
959  }
960  */
961 
962 #ifdef DEBUG
963  if (debug && g_debug) {
964  for (int n = 0; n < N_proc; ++n) {
965  dmutex_guard;
966  std::cout << n << ": jacobian" << std::endl;
967  printf("%5f %5f %5f %5f %5f %5f\n",
968  errorProp(n, 0, 0),
969  errorProp(n, 0, 1),
970  errorProp(n, 0, 2),
971  errorProp(n, 0, 3),
972  errorProp(n, 0, 4),
973  errorProp(n, 0, 5));
974  printf("%5f %5f %5f %5f %5f %5f\n",
975  errorProp(n, 1, 0),
976  errorProp(n, 1, 1),
977  errorProp(n, 1, 2),
978  errorProp(n, 1, 3),
979  errorProp(n, 1, 4),
980  errorProp(n, 1, 5));
981  printf("%5f %5f %5f %5f %5f %5f\n",
982  errorProp(n, 2, 0),
983  errorProp(n, 2, 1),
984  errorProp(n, 2, 2),
985  errorProp(n, 2, 3),
986  errorProp(n, 2, 4),
987  errorProp(n, 2, 5));
988  printf("%5f %5f %5f %5f %5f %5f\n",
989  errorProp(n, 3, 0),
990  errorProp(n, 3, 1),
991  errorProp(n, 3, 2),
992  errorProp(n, 3, 3),
993  errorProp(n, 3, 4),
994  errorProp(n, 3, 5));
995  printf("%5f %5f %5f %5f %5f %5f\n",
996  errorProp(n, 4, 0),
997  errorProp(n, 4, 1),
998  errorProp(n, 4, 2),
999  errorProp(n, 4, 3),
1000  errorProp(n, 4, 4),
1001  errorProp(n, 4, 5));
1002  printf("%5f %5f %5f %5f %5f %5f\n",
1003  errorProp(n, 5, 0),
1004  errorProp(n, 5, 1),
1005  errorProp(n, 5, 2),
1006  errorProp(n, 5, 3),
1007  errorProp(n, 5, 4),
1008  errorProp(n, 5, 5));
1009  }
1010  }
1011 #endif
1012  }
1013 
1014  //==============================================================================
1015 
1016  void applyMaterialEffects(const MPlexQF& hitsRl,
1017  const MPlexQF& hitsXi,
1018  const MPlexQF& propSign,
1019  MPlexLS& outErr,
1020  MPlexLV& outPar,
1021  const int N_proc,
1022  const bool isBarrel) {
1023 #pragma omp simd
1024  for (int n = 0; n < NN; ++n) {
1025  float radL = hitsRl.constAt(n, 0, 0);
1026  if (radL < 1e-13f)
1027  continue; //ugly, please fixme
1028  const float theta = outPar.constAt(n, 5, 0);
1029  const float pt = 1.f / outPar.constAt(n, 3, 0); //fixme, make sure it is positive?
1030  const float p = pt / std::sin(theta);
1031  const float p2 = p * p;
1032  constexpr float mpi = 0.140; // m=140 MeV, pion
1033  constexpr float mpi2 = mpi * mpi; // m=140 MeV, pion
1034  const float beta2 = p2 / (p2 + mpi2);
1035  const float beta = std::sqrt(beta2);
1036  //radiation lenght, corrected for the crossing angle (cos alpha from dot product of radius vector and momentum)
1037  const float invCos = (isBarrel ? p / pt : 1.f / std::abs(std::cos(theta)));
1038  radL = radL * invCos; //fixme works only for barrel geom
1039  // multiple scattering
1040  //vary independently phi and theta by the rms of the planar multiple scattering angle
1041  // XXX-KMD radL < 0, see your fixme above! Repeating bailout
1042  if (radL < 1e-13f)
1043  continue;
1044  // const float thetaMSC = 0.0136f*std::sqrt(radL)*(1.f+0.038f*std::log(radL))/(beta*p);// eq 32.15
1045  // const float thetaMSC2 = thetaMSC*thetaMSC;
1046  const float thetaMSC = 0.0136f * (1.f + 0.038f * std::log(radL)) / (beta * p); // eq 32.15
1047  const float thetaMSC2 = thetaMSC * thetaMSC * radL;
1048  outErr.At(n, 4, 4) += thetaMSC2;
1049  // outErr.At(n, 4, 5) += thetaMSC2;
1050  outErr.At(n, 5, 5) += thetaMSC2;
1051  //std::cout << "beta=" << beta << " p=" << p << std::endl;
1052  //std::cout << "multiple scattering thetaMSC=" << thetaMSC << " thetaMSC2=" << thetaMSC2 << " radL=" << radL << std::endl;
1053  // energy loss
1054  // XXX-KMD beta2 = 1 => 1 / sqrt(0)
1055  // const float gamma = 1.f/std::sqrt(1.f - std::min(beta2, 0.999999f));
1056  // const float gamma2 = gamma*gamma;
1057  const float gamma2 = (p2 + mpi2) / mpi2;
1058  const float gamma = std::sqrt(gamma2); //1.f/std::sqrt(1.f - std::min(beta2, 0.999999f));
1059  constexpr float me = 0.0005; // m=0.5 MeV, electron
1060  const float wmax = 2.f * me * beta2 * gamma2 / (1.f + 2.f * gamma * me / mpi + me * me / (mpi * mpi));
1061  constexpr float I = 16.0e-9 * 10.75;
1062  const float deltahalf = std::log(28.816e-9f * std::sqrt(2.33f * 0.498f) / I) + std::log(beta * gamma) - 0.5f;
1063  const float dEdx =
1064  beta < 1.f
1065  ? (2.f * (hitsXi.constAt(n, 0, 0) * invCos *
1066  (0.5f * std::log(2.f * me * beta2 * gamma2 * wmax / (I * I)) - beta2 - deltahalf) / beta2))
1067  : 0.f; //protect against infs and nans
1068  // dEdx = dEdx*2.;//xi in cmssw is defined with an extra factor 0.5 with respect to formula 27.1 in pdg
1069  //std::cout << "dEdx=" << dEdx << " delta=" << deltahalf << " wmax=" << wmax << " Xi=" << hitsXi.constAt(n,0,0) << std::endl;
1070  const float dP = propSign.constAt(n, 0, 0) * dEdx / beta;
1071  outPar.At(n, 3, 0) = p / (std::max(p + dP, 0.001f) * pt); //stay above 1MeV
1072  //assume 100% uncertainty
1073  outErr.At(n, 3, 3) += dP * dP / (p2 * pt * pt);
1074  }
1075  }
1076 
1077 } // 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:429
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:96
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:594
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
l1tTrackerHTMiss_cfi.deltaZ
deltaZ
Definition: l1tTrackerHTMiss_cfi.py:13
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:487
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
HltBtagPostValidation_cff.c
c
Definition: HltBtagPostValidation_cff.py:31
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
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
mkfit::g_debug
bool g_debug
Definition: Debug.cc:2
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:96
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::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:708
mkfit::Config::useTrigApprox
constexpr bool useTrigApprox
Definition: Config.h:93
correctionTermsCaloMet_cff.C
C
Definition: correctionTermsCaloMet_cff.py:34
HLT_2023v12_cff.beta
beta
Definition: HLT_2023v12_cff.py:7832
debug
#define debug
Definition: HDRShower.cc:19
mkfit::Const::sol
constexpr float sol
Definition: Config.h:56
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:167
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:1016
GetRecoTauVFromDQM_MC_cff.kk
kk
Definition: GetRecoTauVFromDQM_MC_cff.py:84
a
double a
Definition: hdecay.h:119
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
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:98
AlCaHLTBitMon_ParallelJobs.p
def p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
mkfit::Config::Niter
constexpr int Niter
Definition: Config.h:92
cms::cuda::wmax
__host__ __device__ V V wmax
Definition: HistoContainer.h:85
Generated for CMSSW Reference Manual by doxygen 1.8.14