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