CMS 3D CMS Logo
Quintuplet.h
Go to the documentation of this file.
1 #ifndef RecoTracker_LSTCore_src_alpaka_Quintuplet_h
2 #define RecoTracker_LSTCore_src_alpaka_Quintuplet_h
3 
4 #include "HeterogeneousCore/AlpakaInterface/interface/workdivision.h"
5 
6 #include "RecoTracker/LSTCore/interface/ObjectRangesSoA.h"
7 #include "RecoTracker/LSTCore/interface/MiniDoubletsSoA.h"
8 #include "RecoTracker/LSTCore/interface/SegmentsSoA.h"
9 #include "RecoTracker/LSTCore/interface/TripletsSoA.h"
10 #include "RecoTracker/LSTCore/interface/QuintupletsSoA.h"
11 #include "RecoTracker/LSTCore/interface/alpaka/Common.h"
12 #include "RecoTracker/LSTCore/interface/ModulesSoA.h"
13 #include "RecoTracker/LSTCore/interface/EndcapGeometry.h"
14 #include "RecoTracker/LSTCore/interface/ObjectRangesSoA.h"
15 
16 #include "NeuralNetwork.h"
17 #include "Hit.h"
18 #include "Triplet.h" // FIXME: need to refactor common functions to a common place
19 
20 namespace ALPAKA_ACCELERATOR_NAMESPACE::lst {
21  ALPAKA_FN_ACC ALPAKA_FN_INLINE bool checkIntervalOverlap(float firstMin,
22  float firstMax,
23  float secondMin,
24  float secondMax) {
25  return ((firstMin <= secondMin) && (secondMin < firstMax)) || ((secondMin < firstMin) && (firstMin < secondMax));
26  }
27 
28  ALPAKA_FN_ACC ALPAKA_FN_INLINE void addQuintupletToMemory(TripletsConst triplets,
29  Quintuplets quintuplets,
30  unsigned int innerTripletIndex,
31  unsigned int outerTripletIndex,
32  uint16_t lowerModule1,
33  uint16_t lowerModule2,
34  uint16_t lowerModule3,
35  uint16_t lowerModule4,
36  uint16_t lowerModule5,
37  float innerRadius,
38  float bridgeRadius,
39  float outerRadius,
40  float regressionG,
41  float regressionF,
42  float regressionRadius,
43  float rzChiSquared,
44  float rPhiChiSquared,
45  float nonAnchorChiSquared,
46  float pt,
47  float eta,
48  float phi,
49  float scores,
50  uint8_t layer,
51  unsigned int quintupletIndex,
52  bool tightCutFlag) {
53  quintuplets.tripletIndices()[quintupletIndex][0] = innerTripletIndex;
54  quintuplets.tripletIndices()[quintupletIndex][1] = outerTripletIndex;
55 
56  quintuplets.lowerModuleIndices()[quintupletIndex][0] = lowerModule1;
57  quintuplets.lowerModuleIndices()[quintupletIndex][1] = lowerModule2;
58  quintuplets.lowerModuleIndices()[quintupletIndex][2] = lowerModule3;
59  quintuplets.lowerModuleIndices()[quintupletIndex][3] = lowerModule4;
60  quintuplets.lowerModuleIndices()[quintupletIndex][4] = lowerModule5;
61  quintuplets.innerRadius()[quintupletIndex] = __F2H(innerRadius);
62  quintuplets.outerRadius()[quintupletIndex] = __F2H(outerRadius);
63  quintuplets.pt()[quintupletIndex] = __F2H(pt);
64  quintuplets.eta()[quintupletIndex] = __F2H(eta);
65  quintuplets.phi()[quintupletIndex] = __F2H(phi);
66  quintuplets.score_rphisum()[quintupletIndex] = __F2H(scores);
67  quintuplets.isDup()[quintupletIndex] = 0;
68  quintuplets.tightCutFlag()[quintupletIndex] = tightCutFlag;
69  quintuplets.regressionRadius()[quintupletIndex] = regressionRadius;
70  quintuplets.regressionG()[quintupletIndex] = regressionG;
71  quintuplets.regressionF()[quintupletIndex] = regressionF;
72  quintuplets.logicalLayers()[quintupletIndex][0] = triplets.logicalLayers()[innerTripletIndex][0];
73  quintuplets.logicalLayers()[quintupletIndex][1] = triplets.logicalLayers()[innerTripletIndex][1];
74  quintuplets.logicalLayers()[quintupletIndex][2] = triplets.logicalLayers()[innerTripletIndex][2];
75  quintuplets.logicalLayers()[quintupletIndex][3] = triplets.logicalLayers()[outerTripletIndex][1];
76  quintuplets.logicalLayers()[quintupletIndex][4] = triplets.logicalLayers()[outerTripletIndex][2];
77 
78  quintuplets.hitIndices()[quintupletIndex][0] = triplets.hitIndices()[innerTripletIndex][0];
79  quintuplets.hitIndices()[quintupletIndex][1] = triplets.hitIndices()[innerTripletIndex][1];
80  quintuplets.hitIndices()[quintupletIndex][2] = triplets.hitIndices()[innerTripletIndex][2];
81  quintuplets.hitIndices()[quintupletIndex][3] = triplets.hitIndices()[innerTripletIndex][3];
82  quintuplets.hitIndices()[quintupletIndex][4] = triplets.hitIndices()[innerTripletIndex][4];
83  quintuplets.hitIndices()[quintupletIndex][5] = triplets.hitIndices()[innerTripletIndex][5];
84  quintuplets.hitIndices()[quintupletIndex][6] = triplets.hitIndices()[outerTripletIndex][2];
85  quintuplets.hitIndices()[quintupletIndex][7] = triplets.hitIndices()[outerTripletIndex][3];
86  quintuplets.hitIndices()[quintupletIndex][8] = triplets.hitIndices()[outerTripletIndex][4];
87  quintuplets.hitIndices()[quintupletIndex][9] = triplets.hitIndices()[outerTripletIndex][5];
88  quintuplets.bridgeRadius()[quintupletIndex] = bridgeRadius;
89  quintuplets.rzChiSquared()[quintupletIndex] = rzChiSquared;
90  quintuplets.chiSquared()[quintupletIndex] = rPhiChiSquared;
91  quintuplets.nonAnchorChiSquared()[quintupletIndex] = nonAnchorChiSquared;
92  }
93 
94  //90% constraint
95  ALPAKA_FN_ACC ALPAKA_FN_INLINE bool passChiSquaredConstraint(ModulesConst modules,
96  uint16_t lowerModuleIndex1,
97  uint16_t lowerModuleIndex2,
98  uint16_t lowerModuleIndex3,
99  uint16_t lowerModuleIndex4,
100  uint16_t lowerModuleIndex5,
101  float chiSquared) {
102  // Using lstLayer numbering convention defined in ModuleMethods.h
103  const int layer1 = modules.lstLayers()[lowerModuleIndex1];
104  const int layer2 = modules.lstLayers()[lowerModuleIndex2];
105  const int layer3 = modules.lstLayers()[lowerModuleIndex3];
106  const int layer4 = modules.lstLayers()[lowerModuleIndex4];
107  const int layer5 = modules.lstLayers()[lowerModuleIndex5];
108 
109  if (layer1 == 7 and layer2 == 8 and layer3 == 9) {
110  if (layer4 == 10 and layer5 == 11) {
111  return chiSquared < 0.01788f;
112  } else if (layer4 == 10 and layer5 == 16) {
113  return chiSquared < 0.04725f;
114  } else if (layer4 == 15 and layer5 == 16) {
115  return chiSquared < 0.04725f;
116  }
117  } else if (layer1 == 1 and layer2 == 7 and layer3 == 8) {
118  if (layer4 == 9 and layer5 == 10) {
119  return chiSquared < 0.01788f;
120  } else if (layer4 == 9 and layer5 == 15) {
121  return chiSquared < 0.08234f;
122  }
123  } else if (layer1 == 1 and layer2 == 2 and layer3 == 7) {
124  if (layer4 == 8 and layer5 == 9) {
125  return chiSquared < 0.02360f;
126  } else if (layer4 == 8 and layer5 == 14) {
127  return chiSquared < 0.07167f;
128  } else if (layer4 == 13 and layer5 == 14) {
129  return chiSquared < 0.08234f;
130  }
131  } else if (layer1 == 1 and layer2 == 2 and layer3 == 3) {
132  if (layer4 == 7 and layer5 == 8) {
133  return chiSquared < 0.01026f;
134  } else if (layer4 == 7 and layer5 == 13) {
135  return chiSquared < 0.06238f;
136  } else if (layer4 == 12 and layer5 == 13) {
137  return chiSquared < 0.06238f;
138  }
139  } else if (layer1 == 1 and layer2 == 2 and layer3 == 3 and layer4 == 4) {
140  if (layer5 == 5) {
141  return chiSquared < 0.04725f;
142  } else if (layer5 == 12) {
143  return chiSquared < 0.09461f;
144  }
145  } else if (layer1 == 2 and layer2 == 7 and layer3 == 8) {
146  if (layer4 == 9 and layer5 == 10) {
147  return chiSquared < 0.00512f;
148  }
149  if (layer4 == 9 and layer5 == 15) {
150  return chiSquared < 0.04112f;
151  } else if (layer4 == 14 and layer5 == 15) {
152  return chiSquared < 0.06238f;
153  }
154  } else if (layer1 == 2 and layer2 == 3 and layer3 == 7) {
155  if (layer4 == 8 and layer5 == 14) {
156  return chiSquared < 0.07167f;
157  } else if (layer4 == 13 and layer5 == 14) {
158  return chiSquared < 0.06238f;
159  }
160  } else if (layer1 == 2 and layer2 == 3 and layer3 == 4) {
161  if (layer4 == 5 and layer5 == 6) {
162  return chiSquared < 0.08234f;
163  } else if (layer4 == 5 and layer5 == 12) {
164  return chiSquared < 0.10870f;
165  } else if (layer4 == 12 and layer5 == 13) {
166  return chiSquared < 0.10870f;
167  }
168  } else if (layer1 == 3 and layer2 == 7 and layer3 == 8 and layer4 == 14 and layer5 == 15) {
169  return chiSquared < 0.09461f;
170  } else if (layer1 == 3 and layer2 == 4 and layer3 == 5 and layer4 == 12 and layer5 == 13) {
171  return chiSquared < 0.09461f;
172  }
173 
174  return true;
175  }
176 
177  //bounds can be found at http://uaf-10.t2.ucsd.edu/~bsathian/SDL/T5_RZFix/t5_rz_thresholds.txt
178  template <typename TAcc>
179  ALPAKA_FN_ACC ALPAKA_FN_INLINE bool passT5RZConstraint(TAcc const& acc,
180  ModulesConst modules,
181  MiniDoubletsConst mds,
182  unsigned int firstMDIndex,
183  unsigned int secondMDIndex,
184  unsigned int thirdMDIndex,
185  unsigned int fourthMDIndex,
186  unsigned int fifthMDIndex,
187  uint16_t lowerModuleIndex1,
188  uint16_t lowerModuleIndex2,
189  uint16_t lowerModuleIndex3,
190  uint16_t lowerModuleIndex4,
191  uint16_t lowerModuleIndex5,
192  float& rzChiSquared,
193  float inner_pt,
194  float innerRadius,
195  float g,
196  float f,
197  bool& tightCutFlag) {
198  //(g,f) is the center of the circle fitted by the innermost 3 points on x,y coordinates
199  const float rt1 = mds.anchorRt()[firstMDIndex] / 100; //in the unit of m instead of cm
200  const float rt2 = mds.anchorRt()[secondMDIndex] / 100;
201  const float rt3 = mds.anchorRt()[thirdMDIndex] / 100;
202  const float rt4 = mds.anchorRt()[fourthMDIndex] / 100;
203  const float rt5 = mds.anchorRt()[fifthMDIndex] / 100;
204 
205  const float z1 = mds.anchorZ()[firstMDIndex] / 100;
206  const float z2 = mds.anchorZ()[secondMDIndex] / 100;
207  const float z3 = mds.anchorZ()[thirdMDIndex] / 100;
208  const float z4 = mds.anchorZ()[fourthMDIndex] / 100;
209  const float z5 = mds.anchorZ()[fifthMDIndex] / 100;
210 
211  // Using lst_layer numbering convention defined in ModuleMethods.h
212  const int layer1 = modules.lstLayers()[lowerModuleIndex1];
213  const int layer2 = modules.lstLayers()[lowerModuleIndex2];
214  const int layer3 = modules.lstLayers()[lowerModuleIndex3];
215  const int layer4 = modules.lstLayers()[lowerModuleIndex4];
216  const int layer5 = modules.lstLayers()[lowerModuleIndex5];
217 
218  //slope computed using the internal T3s
219  const int moduleType1 = modules.moduleType()[lowerModuleIndex1]; //0 is ps, 1 is 2s
220  const int moduleType2 = modules.moduleType()[lowerModuleIndex2];
221  const int moduleType3 = modules.moduleType()[lowerModuleIndex3];
222  const int moduleType4 = modules.moduleType()[lowerModuleIndex4];
223  const int moduleType5 = modules.moduleType()[lowerModuleIndex5];
224 
225  const float x1 = mds.anchorX()[firstMDIndex] / 100;
226  const float x2 = mds.anchorX()[secondMDIndex] / 100;
227  const float x3 = mds.anchorX()[thirdMDIndex] / 100;
228  const float x4 = mds.anchorX()[fourthMDIndex] / 100;
229  const float y1 = mds.anchorY()[firstMDIndex] / 100;
230  const float y2 = mds.anchorY()[secondMDIndex] / 100;
231  const float y3 = mds.anchorY()[thirdMDIndex] / 100;
232  const float y4 = mds.anchorY()[fourthMDIndex] / 100;
233 
234  float residual = 0;
235  float error2 = 0;
236  float x_center = g / 100, y_center = f / 100;
237  float x_init = mds.anchorX()[thirdMDIndex] / 100;
238  float y_init = mds.anchorY()[thirdMDIndex] / 100;
239  float z_init = mds.anchorZ()[thirdMDIndex] / 100;
240  float rt_init = mds.anchorRt()[thirdMDIndex] / 100; //use the second MD as initial point
241 
242  if (moduleType3 == 1) // 1: if MD3 is in 2s layer
243  {
244  x_init = mds.anchorX()[secondMDIndex] / 100;
245  y_init = mds.anchorY()[secondMDIndex] / 100;
246  z_init = mds.anchorZ()[secondMDIndex] / 100;
247  rt_init = mds.anchorRt()[secondMDIndex] / 100;
248  }
249 
250  // start from a circle of inner T3.
251  // to determine the charge
252  int charge = 0;
253  float slope3c = (y3 - y_center) / (x3 - x_center);
254  float slope1c = (y1 - y_center) / (x1 - x_center);
255  // these 4 "if"s basically separate the x-y plane into 4 quarters. It determines geometrically how a circle and line slope goes and their positions, and we can get the charges correspondingly.
256  if ((y3 - y_center) > 0 && (y1 - y_center) > 0) {
257  if (slope1c > 0 && slope3c < 0)
258  charge = -1; // on x axis of a quarter, 3 hits go anti-clockwise
259  else if (slope1c < 0 && slope3c > 0)
260  charge = 1; // on x axis of a quarter, 3 hits go clockwise
261  else if (slope3c > slope1c)
262  charge = -1;
263  else if (slope3c < slope1c)
264  charge = 1;
265  } else if ((y3 - y_center) < 0 && (y1 - y_center) < 0) {
266  if (slope1c < 0 && slope3c > 0)
267  charge = 1;
268  else if (slope1c > 0 && slope3c < 0)
269  charge = -1;
270  else if (slope3c > slope1c)
271  charge = -1;
272  else if (slope3c < slope1c)
273  charge = 1;
274  } else if ((y3 - y_center) < 0 && (y1 - y_center) > 0) {
275  if ((x3 - x_center) > 0 && (x1 - x_center) > 0)
276  charge = 1;
277  else if ((x3 - x_center) < 0 && (x1 - x_center) < 0)
278  charge = -1;
279  } else if ((y3 - y_center) > 0 && (y1 - y_center) < 0) {
280  if ((x3 - x_center) > 0 && (x1 - x_center) > 0)
281  charge = -1;
282  else if ((x3 - x_center) < 0 && (x1 - x_center) < 0)
283  charge = 1;
284  }
285 
286  float pseudo_phi = alpaka::math::atan(
287  acc, (y_init - y_center) / (x_init - x_center)); //actually represent pi/2-phi, wrt helix axis z
288  float Pt = inner_pt, Px = Pt * alpaka::math::abs(acc, alpaka::math::sin(acc, pseudo_phi)),
289  Py = Pt * alpaka::math::abs(acc, cos(pseudo_phi));
290 
291  // Above line only gives you the correct value of Px and Py, but signs of Px and Py calculated below.
292  // We look at if the circle is clockwise or anti-clock wise, to make it simpler, we separate the x-y plane into 4 quarters.
293  if (x_init > x_center && y_init > y_center) //1st quad
294  {
295  if (charge == 1)
296  Py = -Py;
297  if (charge == -1)
298  Px = -Px;
299  }
300  if (x_init < x_center && y_init > y_center) //2nd quad
301  {
302  if (charge == -1) {
303  Px = -Px;
304  Py = -Py;
305  }
306  }
307  if (x_init < x_center && y_init < y_center) //3rd quad
308  {
309  if (charge == 1)
310  Px = -Px;
311  if (charge == -1)
312  Py = -Py;
313  }
314  if (x_init > x_center && y_init < y_center) //4th quad
315  {
316  if (charge == 1) {
317  Px = -Px;
318  Py = -Py;
319  }
320  }
321 
322  // But if the initial T5 curve goes across quarters(i.e. cross axis to separate the quarters), need special redeclaration of Px,Py signs on these to avoid errors
323  if (moduleType3 == 0) { // 0 is ps
324  if (x4 < x3 && x3 < x2)
325  Px = -alpaka::math::abs(acc, Px);
326  else if (x4 > x3 && x3 > x2)
327  Px = alpaka::math::abs(acc, Px);
328  if (y4 < y3 && y3 < y2)
329  Py = -alpaka::math::abs(acc, Py);
330  else if (y4 > y3 && y3 > y2)
331  Py = alpaka::math::abs(acc, Py);
332  } else if (moduleType3 == 1) // 1 is 2s
333  {
334  if (x3 < x2 && x2 < x1)
335  Px = -alpaka::math::abs(acc, Px);
336  else if (x3 > x2 && x2 > x1)
337  Px = alpaka::math::abs(acc, Px);
338  if (y3 < y2 && y2 < y1)
339  Py = -alpaka::math::abs(acc, Py);
340  else if (y3 > y2 && y2 > y1)
341  Py = alpaka::math::abs(acc, Py);
342  }
343 
344  //to get Pz, we use pt/pz=ds/dz, ds is the arclength between MD1 and MD3.
345  float AO = alpaka::math::sqrt(acc, (x1 - x_center) * (x1 - x_center) + (y1 - y_center) * (y1 - y_center));
346  float BO =
347  alpaka::math::sqrt(acc, (x_init - x_center) * (x_init - x_center) + (y_init - y_center) * (y_init - y_center));
348  float AB2 = (x1 - x_init) * (x1 - x_init) + (y1 - y_init) * (y1 - y_init);
349  float dPhi = alpaka::math::acos(acc, (AO * AO + BO * BO - AB2) / (2 * AO * BO));
350  float ds = innerRadius / 100 * dPhi;
351 
352  float Pz = (z_init - z1) / ds * Pt;
353  float p = alpaka::math::sqrt(acc, Px * Px + Py * Py + Pz * Pz);
354 
355  float a = -2.f * k2Rinv1GeVf * 100 * charge; // multiply by 100 to make the correct length units
356 
357  float zsi, rtsi;
358  int layeri, moduleTypei;
359  rzChiSquared = 0;
360  for (size_t i = 2; i < 6; i++) {
361  if (i == 2) {
362  zsi = z2;
363  rtsi = rt2;
364  layeri = layer2;
365  moduleTypei = moduleType2;
366  } else if (i == 3) {
367  zsi = z3;
368  rtsi = rt3;
369  layeri = layer3;
370  moduleTypei = moduleType3;
371  } else if (i == 4) {
372  zsi = z4;
373  rtsi = rt4;
374  layeri = layer4;
375  moduleTypei = moduleType4;
376  } else if (i == 5) {
377  zsi = z5;
378  rtsi = rt5;
379  layeri = layer5;
380  moduleTypei = moduleType5;
381  }
382 
383  if (moduleType3 == 0) { //0: ps
384  if (i == 3)
385  continue;
386  } else {
387  if (i == 2)
388  continue;
389  }
390 
391  // calculation is copied from PixelTriplet.h computePT3RZChiSquared
392  float diffr = 0, diffz = 0;
393 
394  float rou = a / p;
395  // for endcap
396  float s = (zsi - z_init) * p / Pz;
397  float x = x_init + Px / a * alpaka::math::sin(acc, rou * s) - Py / a * (1 - alpaka::math::cos(acc, rou * s));
398  float y = y_init + Py / a * alpaka::math::sin(acc, rou * s) + Px / a * (1 - alpaka::math::cos(acc, rou * s));
399  diffr = (rtsi - alpaka::math::sqrt(acc, x * x + y * y)) * 100;
400 
401  // for barrel
402  if (layeri <= 6) {
403  float paraA =
404  rt_init * rt_init + 2 * (Px * Px + Py * Py) / (a * a) + 2 * (y_init * Px - x_init * Py) / a - rtsi * rtsi;
405  float paraB = 2 * (x_init * Px + y_init * Py) / a;
406  float paraC = 2 * (y_init * Px - x_init * Py) / a + 2 * (Px * Px + Py * Py) / (a * a);
407  float A = paraB * paraB + paraC * paraC;
408  float B = 2 * paraA * paraB;
409  float C = paraA * paraA - paraC * paraC;
410  float sol1 = (-B + alpaka::math::sqrt(acc, B * B - 4 * A * C)) / (2 * A);
411  float sol2 = (-B - alpaka::math::sqrt(acc, B * B - 4 * A * C)) / (2 * A);
412  float solz1 = alpaka::math::asin(acc, sol1) / rou * Pz / p + z_init;
413  float solz2 = alpaka::math::asin(acc, sol2) / rou * Pz / p + z_init;
414  float diffz1 = (solz1 - zsi) * 100;
415  float diffz2 = (solz2 - zsi) * 100;
416  if (alpaka::math::isnan(acc, diffz1))
417  diffz = diffz2;
418  else if (alpaka::math::isnan(acc, diffz2))
419  diffz = diffz1;
420  else {
421  diffz = (alpaka::math::abs(acc, diffz1) < alpaka::math::abs(acc, diffz2)) ? diffz1 : diffz2;
422  }
423  }
424  residual = (layeri > 6) ? diffr : diffz;
425 
426  //PS Modules
427  if (moduleTypei == 0) {
428  error2 = kPixelPSZpitch * kPixelPSZpitch;
429  } else //2S modules
430  {
431  error2 = kStrip2SZpitch * kStrip2SZpitch;
432  }
433 
434  //check the tilted module, side: PosZ, NegZ, Center(for not tilted)
435  float drdz;
436  short side, subdets;
437  if (i == 2) {
438  drdz = alpaka::math::abs(acc, modules.drdzs()[lowerModuleIndex2]);
439  side = modules.sides()[lowerModuleIndex2];
440  subdets = modules.subdets()[lowerModuleIndex2];
441  }
442  if (i == 3) {
443  drdz = alpaka::math::abs(acc, modules.drdzs()[lowerModuleIndex3]);
444  side = modules.sides()[lowerModuleIndex3];
445  subdets = modules.subdets()[lowerModuleIndex3];
446  }
447  if (i == 2 || i == 3) {
448  residual = (layeri <= 6 && ((side == Center) or (drdz < 1))) ? diffz : diffr;
449  float projection_missing2 = 1.f;
450  if (drdz < 1)
451  projection_missing2 =
452  ((subdets == Endcap) or (side == Center)) ? 1.f : 1.f / (1 + drdz * drdz); // cos(atan(drdz)), if dr/dz<1
453  if (drdz > 1)
454  projection_missing2 = ((subdets == Endcap) or (side == Center))
455  ? 1.f
456  : (drdz * drdz) / (1 + drdz * drdz); //sin(atan(drdz)), if dr/dz>1
457  error2 = error2 * projection_missing2;
458  }
459  rzChiSquared += 12 * (residual * residual) / error2;
460  }
461  // for set rzchi2 cut
462  // if the 5 points are linear, helix calculation gives nan
463  if (inner_pt > 100 || alpaka::math::isnan(acc, rzChiSquared)) {
464  float slope;
465  if (moduleType1 == 0 and moduleType2 == 0 and moduleType3 == 1) //PSPS2S
466  {
467  slope = (z2 - z1) / (rt2 - rt1);
468  } else {
469  slope = (z3 - z1) / (rt3 - rt1);
470  }
471  float residual4_linear = (layer4 <= 6) ? ((z4 - z1) - slope * (rt4 - rt1)) : ((rt4 - rt1) - (z4 - z1) / slope);
472  float residual5_linear = (layer4 <= 6) ? ((z5 - z1) - slope * (rt5 - rt1)) : ((rt5 - rt1) - (z5 - z1) / slope);
473 
474  // creating a chi squared type quantity
475  // 0-> PS, 1->2S
476  residual4_linear = (moduleType4 == 0) ? residual4_linear / kPixelPSZpitch : residual4_linear / kStrip2SZpitch;
477  residual5_linear = (moduleType5 == 0) ? residual5_linear / kPixelPSZpitch : residual5_linear / kStrip2SZpitch;
478  residual4_linear = residual4_linear * 100;
479  residual5_linear = residual5_linear * 100;
480 
481  rzChiSquared = 12 * (residual4_linear * residual4_linear + residual5_linear * residual5_linear);
482  return rzChiSquared < 4.677f;
483  }
484 
485  // when building T5, apply 99% chi2 cuts as default, and add to pT5 collection. But when adding T5 to TC collections, apply 95% cut to reduce the fake rate
486  tightCutFlag = false;
487  // The category numbers are related to module regions and layers, decoding of the region numbers can be found here in slide 2 table. https://github.com/SegmentLinking/TrackLooper/files/11420927/part.2.pdf
488  // The commented numbers after each case is the region code, and can look it up from the table to see which category it belongs to. For example, //0 means T5 built with Endcap 1,2,3,4,5 ps modules
489  if (layer1 == 7 and layer2 == 8 and layer3 == 9 and layer4 == 10 and layer5 == 11) //0
490  {
491  if (rzChiSquared < 94.470f)
492  tightCutFlag = true;
493  return true;
494  } else if (layer1 == 7 and layer2 == 8 and layer3 == 9 and layer4 == 10 and layer5 == 16) //1
495  {
496  if (rzChiSquared < 22.099f)
497  tightCutFlag = true;
498  return rzChiSquared < 37.956f;
499  } else if (layer1 == 7 and layer2 == 8 and layer3 == 9 and layer4 == 15 and layer5 == 16) //2
500  {
501  if (rzChiSquared < 7.992f)
502  tightCutFlag = true;
503  return rzChiSquared < 11.622f;
504  } else if (layer1 == 1 and layer2 == 7 and layer3 == 8 and layer4 == 9) {
505  if (layer5 == 10) //3
506  {
507  if (rzChiSquared < 111.390f)
508  tightCutFlag = true;
509  return true;
510  }
511  if (layer5 == 15) //4
512  {
513  if (rzChiSquared < 18.351f)
514  tightCutFlag = true;
515  return rzChiSquared < 37.941f;
516  }
517  } else if (layer1 == 1 and layer2 == 2 and layer3 == 7) {
518  if (layer4 == 8 and layer5 == 9) //5
519  {
520  if (rzChiSquared < 116.148f)
521  tightCutFlag = true;
522  return true;
523  }
524  if (layer4 == 8 and layer5 == 14) //6
525  {
526  if (rzChiSquared < 19.352f)
527  tightCutFlag = true;
528  return rzChiSquared < 52.561f;
529  } else if (layer4 == 13 and layer5 == 14) //7
530  {
531  if (rzChiSquared < 10.392f)
532  tightCutFlag = true;
533  return rzChiSquared < 13.76f;
534  }
535  } else if (layer1 == 1 and layer2 == 2 and layer3 == 3) {
536  if (layer4 == 7 and layer5 == 8) //8
537  {
538  if (rzChiSquared < 27.824f)
539  tightCutFlag = true;
540  return rzChiSquared < 44.247f;
541  } else if (layer4 == 7 and layer5 == 13) //9
542  {
543  if (rzChiSquared < 18.145f)
544  tightCutFlag = true;
545  return rzChiSquared < 33.752f;
546  } else if (layer4 == 12 and layer5 == 13) //10
547  {
548  if (rzChiSquared < 13.308f)
549  tightCutFlag = true;
550  return rzChiSquared < 21.213f;
551  } else if (layer4 == 4 and layer5 == 5) //11
552  {
553  if (rzChiSquared < 15.627f)
554  tightCutFlag = true;
555  return rzChiSquared < 29.035f;
556  } else if (layer4 == 4 and layer5 == 12) //12
557  {
558  if (rzChiSquared < 14.64f)
559  tightCutFlag = true;
560  return rzChiSquared < 23.037f;
561  }
562  } else if (layer1 == 2 and layer2 == 7 and layer3 == 8) {
563  if (layer4 == 9 and layer5 == 15) //14
564  {
565  if (rzChiSquared < 24.662f)
566  tightCutFlag = true;
567  return rzChiSquared < 41.036f;
568  } else if (layer4 == 14 and layer5 == 15) //15
569  {
570  if (rzChiSquared < 8.866f)
571  tightCutFlag = true;
572  return rzChiSquared < 14.092f;
573  }
574  } else if (layer1 == 2 and layer2 == 3 and layer3 == 7) {
575  if (layer4 == 8 and layer5 == 14) //16
576  {
577  if (rzChiSquared < 23.730f)
578  tightCutFlag = true;
579  return rzChiSquared < 23.748f;
580  }
581  if (layer4 == 13 and layer5 == 14) //17
582  {
583  if (rzChiSquared < 10.772f)
584  tightCutFlag = true;
585  return rzChiSquared < 17.945f;
586  }
587  } else if (layer1 == 2 and layer2 == 3 and layer3 == 4) {
588  if (layer4 == 5 and layer5 == 6) //18
589  {
590  if (rzChiSquared < 6.065f)
591  tightCutFlag = true;
592  return rzChiSquared < 8.803f;
593  } else if (layer4 == 5 and layer5 == 12) //19
594  {
595  if (rzChiSquared < 5.693f)
596  tightCutFlag = true;
597  return rzChiSquared < 7.930f;
598  }
599 
600  else if (layer4 == 12 and layer5 == 13) //20
601  {
602  if (rzChiSquared < 5.473f)
603  tightCutFlag = true;
604  return rzChiSquared < 7.626f;
605  }
606  }
607  return true;
608  }
609 
610  template <typename TAcc>
611  ALPAKA_FN_ACC ALPAKA_FN_INLINE bool T5HasCommonMiniDoublet(TripletsConst triplets,
612  SegmentsConst segments,
613  unsigned int innerTripletIndex,
614  unsigned int outerTripletIndex) {
615  unsigned int innerOuterSegmentIndex = triplets.segmentIndices()[innerTripletIndex][1];
616  unsigned int outerInnerSegmentIndex = triplets.segmentIndices()[outerTripletIndex][0];
617  unsigned int innerOuterOuterMiniDoubletIndex =
618  segments.mdIndices()[innerOuterSegmentIndex][1]; //inner triplet outer segment outer MD index
619  unsigned int outerInnerInnerMiniDoubletIndex =
620  segments.mdIndices()[outerInnerSegmentIndex][0]; //outer triplet inner segment inner MD index
621 
622  return (innerOuterOuterMiniDoubletIndex == outerInnerInnerMiniDoubletIndex);
623  }
624 
625  template <typename TAcc>
626  ALPAKA_FN_ACC ALPAKA_FN_INLINE void computeErrorInRadius(TAcc const& acc,
627  float* x1Vec,
628  float* y1Vec,
629  float* x2Vec,
630  float* y2Vec,
631  float* x3Vec,
632  float* y3Vec,
633  float& minimumRadius,
634  float& maximumRadius) {
635  //brute force
636  float candidateRadius;
637  float g, f;
638  minimumRadius = kVerticalModuleSlope;
639  maximumRadius = 0.f;
640  for (size_t i = 0; i < 3; i++) {
641  float x1 = x1Vec[i];
642  float y1 = y1Vec[i];
643  for (size_t j = 0; j < 3; j++) {
644  float x2 = x2Vec[j];
645  float y2 = y2Vec[j];
646  for (size_t k = 0; k < 3; k++) {
647  float x3 = x3Vec[k];
648  float y3 = y3Vec[k];
649  candidateRadius = computeRadiusFromThreeAnchorHits(acc, x1, y1, x2, y2, x3, y3, g, f);
650  maximumRadius = alpaka::math::max(acc, candidateRadius, maximumRadius);
651  minimumRadius = alpaka::math::min(acc, candidateRadius, minimumRadius);
652  }
653  }
654  }
655  }
656 
657  template <typename TAcc>
658  ALPAKA_FN_ACC ALPAKA_FN_INLINE bool matchRadiiBBBEE12378(TAcc const& acc,
659  float innerRadius,
660  float bridgeRadius,
661  float outerRadius,
662  float bridgeRadiusMin2S,
663  float bridgeRadiusMax2S) {
664  float innerInvRadiusMin, innerInvRadiusMax, bridgeInvRadiusMin, bridgeInvRadiusMax;
665 
666  float innerInvRadiusErrorBound = 0.178f;
667  float bridgeInvRadiusErrorBound = 0.507f;
668 
669  innerInvRadiusMax = (1.f + innerInvRadiusErrorBound) / innerRadius;
670  innerInvRadiusMin = alpaka::math::max(acc, 0.f, (1.f - innerInvRadiusErrorBound) / innerRadius);
671 
672  bridgeInvRadiusMax = (1.f + bridgeInvRadiusErrorBound) / bridgeRadius;
673  bridgeInvRadiusMin = alpaka::math::max(acc, 0.f, (1.f - bridgeInvRadiusErrorBound) / bridgeRadius);
674 
675  return checkIntervalOverlap(innerInvRadiusMin,
676  innerInvRadiusMax,
677  alpaka::math::min(acc, bridgeInvRadiusMin, 1.0f / bridgeRadiusMax2S),
678  alpaka::math::max(acc, bridgeInvRadiusMax, 1.0f / bridgeRadiusMin2S));
679  }
680 
681  /*bounds for high Pt taken from : http://uaf-10.t2.ucsd.edu/~bsathian/SDL/T5_efficiency/efficiencies/new_efficiencies/efficiencies_20210513_T5_recovering_high_Pt_efficiencies/highE_radius_matching/highE_bounds.txt */
682  template <typename TAcc>
683  ALPAKA_FN_ACC ALPAKA_FN_INLINE bool matchRadiiBBBBB(TAcc const& acc,
684  float innerRadius,
685  float bridgeRadius,
686  float outerRadius) {
687  float innerInvRadiusMin, innerInvRadiusMax, bridgeInvRadiusMin, bridgeInvRadiusMax;
688 
689  float innerInvRadiusErrorBound = 0.1512f;
690  float bridgeInvRadiusErrorBound = 0.1781f;
691 
692  if (innerRadius * k2Rinv1GeVf > 1.f) {
693  innerInvRadiusErrorBound = 0.4449f;
694  bridgeInvRadiusErrorBound = 0.4033f;
695  }
696 
697  innerInvRadiusMax = (1.f + innerInvRadiusErrorBound) / innerRadius;
698  innerInvRadiusMin = alpaka::math::max(acc, 0.f, (1.f - innerInvRadiusErrorBound) / innerRadius);
699 
700  bridgeInvRadiusMax = (1.f + bridgeInvRadiusErrorBound) / bridgeRadius;
701  bridgeInvRadiusMin = alpaka::math::max(acc, 0.f, (1.f - bridgeInvRadiusErrorBound) / bridgeRadius);
702 
703  return checkIntervalOverlap(innerInvRadiusMin, innerInvRadiusMax, bridgeInvRadiusMin, bridgeInvRadiusMax);
704  }
705 
706  template <typename TAcc>
707  ALPAKA_FN_ACC ALPAKA_FN_INLINE bool matchRadiiBBBBE(TAcc const& acc,
708  float innerRadius,
709  float bridgeRadius,
710  float outerRadius) {
711  float innerInvRadiusMin, innerInvRadiusMax, bridgeInvRadiusMin, bridgeInvRadiusMax;
712 
713  float innerInvRadiusErrorBound = 0.1781f;
714  float bridgeInvRadiusErrorBound = 0.2167f;
715 
716  if (innerRadius * k2Rinv1GeVf > 1.f) {
717  innerInvRadiusErrorBound = 0.4750f;
718  bridgeInvRadiusErrorBound = 0.3903f;
719  }
720 
721  innerInvRadiusMax = (1.f + innerInvRadiusErrorBound) / innerRadius;
722  innerInvRadiusMin = alpaka::math::max(acc, 0.f, (1.f - innerInvRadiusErrorBound) / innerRadius);
723 
724  bridgeInvRadiusMax = (1.f + bridgeInvRadiusErrorBound) / bridgeRadius;
725  bridgeInvRadiusMin = alpaka::math::max(acc, 0.f, (1.f - bridgeInvRadiusErrorBound) / bridgeRadius);
726 
727  return checkIntervalOverlap(innerInvRadiusMin, innerInvRadiusMax, bridgeInvRadiusMin, bridgeInvRadiusMax);
728  }
729 
730  template <typename TAcc>
731  ALPAKA_FN_ACC ALPAKA_FN_INLINE bool matchRadiiBBBEE23478(TAcc const& acc,
732  float innerRadius,
733  float bridgeRadius,
734  float outerRadius,
735  float bridgeRadiusMin2S,
736  float bridgeRadiusMax2S) {
737  float innerInvRadiusMin, innerInvRadiusMax, bridgeInvRadiusMin, bridgeInvRadiusMax;
738 
739  float innerInvRadiusErrorBound = 0.2097f;
740  float bridgeInvRadiusErrorBound = 0.8557f;
741 
742  innerInvRadiusMax = (1.f + innerInvRadiusErrorBound) / innerRadius;
743  innerInvRadiusMin = alpaka::math::max(acc, 0.f, (1.f - innerInvRadiusErrorBound) / innerRadius);
744 
745  bridgeInvRadiusMax = (1.f + bridgeInvRadiusErrorBound) / bridgeRadius;
746  bridgeInvRadiusMin = alpaka::math::max(acc, 0.f, (1.f - bridgeInvRadiusErrorBound) / bridgeRadius);
747 
748  return checkIntervalOverlap(innerInvRadiusMin,
749  innerInvRadiusMax,
750  alpaka::math::min(acc, bridgeInvRadiusMin, 1.0f / bridgeRadiusMax2S),
751  alpaka::math::max(acc, bridgeInvRadiusMax, 1.0f / bridgeRadiusMin2S));
752  }
753 
754  template <typename TAcc>
755  ALPAKA_FN_ACC ALPAKA_FN_INLINE bool matchRadiiBBBEE34578(TAcc const& acc,
756  float innerRadius,
757  float bridgeRadius,
758  float outerRadius,
759  float bridgeRadiusMin2S,
760  float bridgeRadiusMax2S) {
761  float innerInvRadiusMin, innerInvRadiusMax, bridgeInvRadiusMin, bridgeInvRadiusMax;
762 
763  float innerInvRadiusErrorBound = 0.066f;
764  float bridgeInvRadiusErrorBound = 0.617f;
765 
766  innerInvRadiusMax = (1.f + innerInvRadiusErrorBound) / innerRadius;
767  innerInvRadiusMin = alpaka::math::max(acc, 0.f, (1.f - innerInvRadiusErrorBound) / innerRadius);
768 
769  bridgeInvRadiusMax = (1.f + bridgeInvRadiusErrorBound) / bridgeRadius;
770  bridgeInvRadiusMin = alpaka::math::max(acc, 0.f, (1.f - bridgeInvRadiusErrorBound) / bridgeRadius);
771 
772  return checkIntervalOverlap(innerInvRadiusMin,
773  innerInvRadiusMax,
774  alpaka::math::min(acc, bridgeInvRadiusMin, 1.0f / bridgeRadiusMax2S),
775  alpaka::math::max(acc, bridgeInvRadiusMax, 1.0f / bridgeRadiusMin2S));
776  }
777 
778  template <typename TAcc>
779  ALPAKA_FN_ACC ALPAKA_FN_INLINE bool matchRadiiBBEEE(TAcc const& acc,
780  float innerRadius,
781  float bridgeRadius,
782  float outerRadius,
783  float bridgeRadiusMin2S,
784  float bridgeRadiusMax2S) {
785  float innerInvRadiusMin, innerInvRadiusMax, bridgeInvRadiusMin, bridgeInvRadiusMax;
786 
787  float innerInvRadiusErrorBound = 0.6376f;
788  float bridgeInvRadiusErrorBound = 2.1381f;
789 
790  if (innerRadius * k2Rinv1GeVf > 1.f) //as good as no selections!
791  {
792  innerInvRadiusErrorBound = 12.9173f;
793  bridgeInvRadiusErrorBound = 5.1700f;
794  }
795 
796  innerInvRadiusMax = (1.f + innerInvRadiusErrorBound) / innerRadius;
797  innerInvRadiusMin = alpaka::math::max(acc, 0.f, (1.f - innerInvRadiusErrorBound) / innerRadius);
798 
799  bridgeInvRadiusMax = (1.f + bridgeInvRadiusErrorBound) / bridgeRadius;
800  bridgeInvRadiusMin = alpaka::math::max(acc, 0.f, (1.f - bridgeInvRadiusErrorBound) / bridgeRadius);
801 
802  return checkIntervalOverlap(innerInvRadiusMin,
803  innerInvRadiusMax,
804  alpaka::math::min(acc, bridgeInvRadiusMin, 1.0f / bridgeRadiusMax2S),
805  alpaka::math::max(acc, bridgeInvRadiusMax, 1.0f / bridgeRadiusMin2S));
806  }
807 
808  template <typename TAcc>
809  ALPAKA_FN_ACC ALPAKA_FN_INLINE bool matchRadiiBEEEE(TAcc const& acc,
810  float innerRadius,
811  float bridgeRadius,
812  float outerRadius,
813  float innerRadiusMin2S,
814  float innerRadiusMax2S,
815  float bridgeRadiusMin2S,
816  float bridgeRadiusMax2S) {
817  float innerInvRadiusMin, innerInvRadiusMax, bridgeInvRadiusMin, bridgeInvRadiusMax;
818 
819  float innerInvRadiusErrorBound = 1.9382f;
820  float bridgeInvRadiusErrorBound = 3.7280f;
821 
822  if (innerRadius * k2Rinv1GeVf > 1.f) {
823  innerInvRadiusErrorBound = 23.2713f;
824  bridgeInvRadiusErrorBound = 21.7980f;
825  }
826 
827  innerInvRadiusMax = (1.f + innerInvRadiusErrorBound) / innerRadius;
828  innerInvRadiusMin = alpaka::math::max(acc, 0.f, (1.f - innerInvRadiusErrorBound) / innerRadius);
829 
830  bridgeInvRadiusMax = (1.f + bridgeInvRadiusErrorBound) / bridgeRadius;
831  bridgeInvRadiusMin = alpaka::math::max(acc, 0.f, (1.f - bridgeInvRadiusErrorBound) / bridgeRadius);
832 
833  return checkIntervalOverlap(alpaka::math::min(acc, innerInvRadiusMin, 1.0f / innerRadiusMax2S),
834  alpaka::math::max(acc, innerInvRadiusMax, 1.0f / innerRadiusMin2S),
835  alpaka::math::min(acc, bridgeInvRadiusMin, 1.0f / bridgeRadiusMax2S),
836  alpaka::math::max(acc, bridgeInvRadiusMax, 1.0f / bridgeRadiusMin2S));
837  }
838 
839  template <typename TAcc>
840  ALPAKA_FN_ACC ALPAKA_FN_INLINE bool matchRadiiEEEEE(TAcc const& acc,
841  float innerRadius,
842  float bridgeRadius,
843  float outerRadius,
844  float innerRadiusMin2S,
845  float innerRadiusMax2S,
846  float bridgeRadiusMin2S,
847  float bridgeRadiusMax2S) {
848  float innerInvRadiusMin, innerInvRadiusMax, bridgeInvRadiusMin, bridgeInvRadiusMax;
849 
850  float innerInvRadiusErrorBound = 1.9382f;
851  float bridgeInvRadiusErrorBound = 2.2091f;
852 
853  if (innerRadius * k2Rinv1GeVf > 1.f) {
854  innerInvRadiusErrorBound = 22.5226f;
855  bridgeInvRadiusErrorBound = 21.0966f;
856  }
857 
858  innerInvRadiusMax = (1.f + innerInvRadiusErrorBound) / innerRadius;
859  innerInvRadiusMin = alpaka::math::max(acc, 0.f, (1.f - innerInvRadiusErrorBound) / innerRadius);
860 
861  bridgeInvRadiusMax = (1.f + bridgeInvRadiusErrorBound) / bridgeRadius;
862  bridgeInvRadiusMin = alpaka::math::max(acc, 0.f, (1.f - bridgeInvRadiusErrorBound) / bridgeRadius);
863 
864  return checkIntervalOverlap(alpaka::math::min(acc, innerInvRadiusMin, 1.0f / innerRadiusMax2S),
865  alpaka::math::max(acc, innerInvRadiusMax, 1.0f / innerRadiusMin2S),
866  alpaka::math::min(acc, bridgeInvRadiusMin, 1.0f / bridgeRadiusMax2S),
867  alpaka::math::max(acc, bridgeInvRadiusMax, 1.0f / bridgeRadiusMin2S));
868  }
869 
870  template <typename TAcc>
871  ALPAKA_FN_ACC ALPAKA_FN_INLINE void computeSigmasForRegression(TAcc const& acc,
872  ModulesConst modules,
873  const uint16_t* lowerModuleIndices,
874  float* delta1,
875  float* delta2,
876  float* slopes,
877  bool* isFlat,
878  unsigned int nPoints = 5,
879  bool anchorHits = true) {
880  /*
881  Bool anchorHits required to deal with a weird edge case wherein
882  the hits ultimately used in the regression are anchor hits, but the
883  lower modules need not all be Pixel Modules (in case of PS). Similarly,
884  when we compute the chi squared for the non-anchor hits, the "partner module"
885  need not always be a PS strip module, but all non-anchor hits sit on strip
886  modules.
887  */
888 
889  ModuleType moduleType;
890  short moduleSubdet, moduleSide;
891  float inv1 = kWidthPS / kWidth2S;
892  float inv2 = kPixelPSZpitch / kWidth2S;
893  float inv3 = kStripPSZpitch / kWidth2S;
894  for (size_t i = 0; i < nPoints; i++) {
895  moduleType = modules.moduleType()[lowerModuleIndices[i]];
896  moduleSubdet = modules.subdets()[lowerModuleIndices[i]];
897  moduleSide = modules.sides()[lowerModuleIndices[i]];
898  const float& drdz = modules.drdzs()[lowerModuleIndices[i]];
899  slopes[i] = modules.dxdys()[lowerModuleIndices[i]];
900  //category 1 - barrel PS flat
901  if (moduleSubdet == Barrel and moduleType == PS and moduleSide == Center) {
902  delta1[i] = inv1;
903  delta2[i] = inv1;
904  slopes[i] = -999.f;
905  isFlat[i] = true;
906  }
907  //category 2 - barrel 2S
908  else if (moduleSubdet == Barrel and moduleType == TwoS) {
909  delta1[i] = 1.f;
910  delta2[i] = 1.f;
911  slopes[i] = -999.f;
912  isFlat[i] = true;
913  }
914  //category 3 - barrel PS tilted
915  else if (moduleSubdet == Barrel and moduleType == PS and moduleSide != Center) {
916  delta1[i] = inv1;
917  isFlat[i] = false;
918 
919  if (anchorHits) {
920  delta2[i] = (inv2 * drdz / alpaka::math::sqrt(acc, 1 + drdz * drdz));
921  } else {
922  delta2[i] = (inv3 * drdz / alpaka::math::sqrt(acc, 1 + drdz * drdz));
923  }
924  }
925  //category 4 - endcap PS
926  else if (moduleSubdet == Endcap and moduleType == PS) {
927  delta1[i] = inv1;
928  isFlat[i] = false;
929 
930  /*
931  despite the type of the module layer of the lower module index,
932  all anchor hits are on the pixel side and all non-anchor hits are
933  on the strip side!
934  */
935  if (anchorHits) {
936  delta2[i] = inv2;
937  } else {
938  delta2[i] = inv3;
939  }
940  }
941  //category 5 - endcap 2S
942  else if (moduleSubdet == Endcap and moduleType == TwoS) {
943  delta1[i] = 1.f;
944  delta2[i] = 500.f * inv1;
945  isFlat[i] = false;
946  } else {
947 #ifdef WARNINGS
948  printf("ERROR!!!!! I SHOULDN'T BE HERE!!!! subdet = %d, type = %d, side = %d\n",
949  moduleSubdet,
950  moduleType,
951  moduleSide);
952 #endif
953  }
954  }
955  }
956 
957  template <typename TAcc>
958  ALPAKA_FN_ACC ALPAKA_FN_INLINE float computeRadiusUsingRegression(TAcc const& acc,
959  unsigned int nPoints,
960  float* xs,
961  float* ys,
962  float* delta1,
963  float* delta2,
964  float* slopes,
965  bool* isFlat,
966  float& g,
967  float& f,
968  float* sigmas2,
969  float& chiSquared) {
970  float radius = 0.f;
971 
972  // Some extra variables
973  // the two variables will be called x1 and x2, and y (which is x^2 + y^2)
974 
975  float sigmaX1Squared = 0.f;
976  float sigmaX2Squared = 0.f;
977  float sigmaX1X2 = 0.f;
978  float sigmaX1y = 0.f;
979  float sigmaX2y = 0.f;
980  float sigmaY = 0.f;
981  float sigmaX1 = 0.f;
982  float sigmaX2 = 0.f;
983  float sigmaOne = 0.f;
984 
985  float xPrime, yPrime, absArctanSlope, angleM;
986  for (size_t i = 0; i < nPoints; i++) {
987  // Computing sigmas is a very tricky affair
988  // if the module is tilted or endcap, we need to use the slopes properly!
989 
990  absArctanSlope = ((slopes[i] != kVerticalModuleSlope) ? alpaka::math::abs(acc, alpaka::math::atan(acc, slopes[i]))
991  : kPi / 2.f);
992 
993  if (xs[i] > 0 and ys[i] > 0) {
994  angleM = kPi / 2.f - absArctanSlope;
995  } else if (xs[i] < 0 and ys[i] > 0) {
996  angleM = absArctanSlope + kPi / 2.f;
997  } else if (xs[i] < 0 and ys[i] < 0) {
998  angleM = -(absArctanSlope + kPi / 2.f);
999  } else if (xs[i] > 0 and ys[i] < 0) {
1000  angleM = -(kPi / 2.f - absArctanSlope);
1001  } else {
1002  angleM = 0;
1003  }
1004 
1005  if (not isFlat[i]) {
1006  xPrime = xs[i] * alpaka::math::cos(acc, angleM) + ys[i] * alpaka::math::sin(acc, angleM);
1007  yPrime = ys[i] * alpaka::math::cos(acc, angleM) - xs[i] * alpaka::math::sin(acc, angleM);
1008  } else {
1009  xPrime = xs[i];
1010  yPrime = ys[i];
1011  }
1012  sigmas2[i] = 4 * ((xPrime * delta1[i]) * (xPrime * delta1[i]) + (yPrime * delta2[i]) * (yPrime * delta2[i]));
1013 
1014  sigmaX1Squared += (xs[i] * xs[i]) / sigmas2[i];
1015  sigmaX2Squared += (ys[i] * ys[i]) / sigmas2[i];
1016  sigmaX1X2 += (xs[i] * ys[i]) / sigmas2[i];
1017  sigmaX1y += (xs[i] * (xs[i] * xs[i] + ys[i] * ys[i])) / sigmas2[i];
1018  sigmaX2y += (ys[i] * (xs[i] * xs[i] + ys[i] * ys[i])) / sigmas2[i];
1019  sigmaY += (xs[i] * xs[i] + ys[i] * ys[i]) / sigmas2[i];
1020  sigmaX1 += xs[i] / sigmas2[i];
1021  sigmaX2 += ys[i] / sigmas2[i];
1022  sigmaOne += 1.0f / sigmas2[i];
1023  }
1024  float denominator = (sigmaX1X2 - sigmaX1 * sigmaX2) * (sigmaX1X2 - sigmaX1 * sigmaX2) -
1025  (sigmaX1Squared - sigmaX1 * sigmaX1) * (sigmaX2Squared - sigmaX2 * sigmaX2);
1026 
1027  float twoG = ((sigmaX2y - sigmaX2 * sigmaY) * (sigmaX1X2 - sigmaX1 * sigmaX2) -
1028  (sigmaX1y - sigmaX1 * sigmaY) * (sigmaX2Squared - sigmaX2 * sigmaX2)) /
1029  denominator;
1030  float twoF = ((sigmaX1y - sigmaX1 * sigmaY) * (sigmaX1X2 - sigmaX1 * sigmaX2) -
1031  (sigmaX2y - sigmaX2 * sigmaY) * (sigmaX1Squared - sigmaX1 * sigmaX1)) /
1032  denominator;
1033 
1034  float c = -(sigmaY - twoG * sigmaX1 - twoF * sigmaX2) / sigmaOne;
1035  g = 0.5f * twoG;
1036  f = 0.5f * twoF;
1037  if (g * g + f * f - c < 0) {
1038 #ifdef WARNINGS
1039  printf("FATAL! r^2 < 0!\n");
1040 #endif
1041  chiSquared = -1;
1042  return -1;
1043  }
1044 
1045  radius = alpaka::math::sqrt(acc, g * g + f * f - c);
1046  // compute chi squared
1047  chiSquared = 0.f;
1048  for (size_t i = 0; i < nPoints; i++) {
1049  chiSquared += (xs[i] * xs[i] + ys[i] * ys[i] - twoG * xs[i] - twoF * ys[i] + c) *
1050  (xs[i] * xs[i] + ys[i] * ys[i] - twoG * xs[i] - twoF * ys[i] + c) / sigmas2[i];
1051  }
1052  return radius;
1053  }
1054 
1055  template <typename TAcc>
1056  ALPAKA_FN_ACC ALPAKA_FN_INLINE float computeChiSquared(TAcc const& acc,
1057  unsigned int nPoints,
1058  float* xs,
1059  float* ys,
1060  float* delta1,
1061  float* delta2,
1062  float* slopes,
1063  bool* isFlat,
1064  float g,
1065  float f,
1066  float radius) {
1067  // given values of (g, f, radius) and a set of points (and its uncertainties)
1068  // compute chi squared
1069  float c = g * g + f * f - radius * radius;
1070  float chiSquared = 0.f;
1071  float absArctanSlope, angleM, xPrime, yPrime, sigma2;
1072  for (size_t i = 0; i < nPoints; i++) {
1073  absArctanSlope = ((slopes[i] != kVerticalModuleSlope) ? alpaka::math::abs(acc, alpaka::math::atan(acc, slopes[i]))
1074  : kPi / 2.f);
1075  if (xs[i] > 0 and ys[i] > 0) {
1076  angleM = kPi / 2.f - absArctanSlope;
1077  } else if (xs[i] < 0 and ys[i] > 0) {
1078  angleM = absArctanSlope + kPi / 2.f;
1079  } else if (xs[i] < 0 and ys[i] < 0) {
1080  angleM = -(absArctanSlope + kPi / 2.f);
1081  } else if (xs[i] > 0 and ys[i] < 0) {
1082  angleM = -(kPi / 2.f - absArctanSlope);
1083  } else {
1084  angleM = 0;
1085  }
1086 
1087  if (not isFlat[i]) {
1088  xPrime = xs[i] * alpaka::math::cos(acc, angleM) + ys[i] * alpaka::math::sin(acc, angleM);
1089  yPrime = ys[i] * alpaka::math::cos(acc, angleM) - xs[i] * alpaka::math::sin(acc, angleM);
1090  } else {
1091  xPrime = xs[i];
1092  yPrime = ys[i];
1093  }
1094  sigma2 = 4 * ((xPrime * delta1[i]) * (xPrime * delta1[i]) + (yPrime * delta2[i]) * (yPrime * delta2[i]));
1095  chiSquared += (xs[i] * xs[i] + ys[i] * ys[i] - 2 * g * xs[i] - 2 * f * ys[i] + c) *
1096  (xs[i] * xs[i] + ys[i] * ys[i] - 2 * g * xs[i] - 2 * f * ys[i] + c) / sigma2;
1097  }
1098  return chiSquared;
1099  }
1100 
1101  template <typename TAcc>
1102  ALPAKA_FN_ACC ALPAKA_FN_INLINE void runDeltaBetaIterationsT5(TAcc const& acc,
1103  float& betaIn,
1104  float& betaOut,
1105  float betaAv,
1106  float& pt_beta,
1107  float sdIn_dr,
1108  float sdOut_dr,
1109  float dr,
1110  float lIn) {
1111  if (lIn == 0) {
1112  betaOut += alpaka::math::copysign(
1113  acc,
1114  alpaka::math::asin(
1115  acc, alpaka::math::min(acc, sdOut_dr * k2Rinv1GeVf / alpaka::math::abs(acc, pt_beta), kSinAlphaMax)),
1116  betaOut);
1117  return;
1118  }
1119 
1120  if (betaIn * betaOut > 0.f and
1121  (alpaka::math::abs(acc, pt_beta) < 4.f * kPt_betaMax or
1122  (lIn >= 11 and alpaka::math::abs(acc, pt_beta) <
1123  8.f * kPt_betaMax))) //and the pt_beta is well-defined; less strict for endcap-endcap
1124  {
1125  const float betaInUpd =
1126  betaIn +
1127  alpaka::math::copysign(
1128  acc,
1129  alpaka::math::asin(
1130  acc, alpaka::math::min(acc, sdIn_dr * k2Rinv1GeVf / alpaka::math::abs(acc, pt_beta), kSinAlphaMax)),
1131  betaIn); //FIXME: need a faster version
1132  const float betaOutUpd =
1133  betaOut +
1134  alpaka::math::copysign(
1135  acc,
1136  alpaka::math::asin(
1137  acc, alpaka::math::min(acc, sdOut_dr * k2Rinv1GeVf / alpaka::math::abs(acc, pt_beta), kSinAlphaMax)),
1138  betaOut); //FIXME: need a faster version
1139  betaAv = 0.5f * (betaInUpd + betaOutUpd);
1140 
1141  //1st update
1142  const float pt_beta_inv =
1143  1.f / alpaka::math::abs(acc, dr * k2Rinv1GeVf / alpaka::math::sin(acc, betaAv)); //get a better pt estimate
1144 
1145  betaIn += alpaka::math::copysign(
1146  acc,
1147  alpaka::math::asin(acc, alpaka::math::min(acc, sdIn_dr * k2Rinv1GeVf * pt_beta_inv, kSinAlphaMax)),
1148  betaIn); //FIXME: need a faster version
1149  betaOut += alpaka::math::copysign(
1150  acc,
1151  alpaka::math::asin(acc, alpaka::math::min(acc, sdOut_dr * k2Rinv1GeVf * pt_beta_inv, kSinAlphaMax)),
1152  betaOut); //FIXME: need a faster version
1153  //update the av and pt
1154  betaAv = 0.5f * (betaIn + betaOut);
1155  //2nd update
1156  pt_beta = dr * k2Rinv1GeVf / alpaka::math::sin(acc, betaAv); //get a better pt estimate
1157  } else if (lIn < 11 && alpaka::math::abs(acc, betaOut) < 0.2f * alpaka::math::abs(acc, betaIn) &&
1158  alpaka::math::abs(acc, pt_beta) < 12.f * kPt_betaMax) //use betaIn sign as ref
1159  {
1160  const float pt_betaIn = dr * k2Rinv1GeVf / alpaka::math::sin(acc, betaIn);
1161 
1162  const float betaInUpd =
1163  betaIn +
1164  alpaka::math::copysign(
1165  acc,
1166  alpaka::math::asin(
1167  acc, alpaka::math::min(acc, sdIn_dr * k2Rinv1GeVf / alpaka::math::abs(acc, pt_betaIn), kSinAlphaMax)),
1168  betaIn); //FIXME: need a faster version
1169  const float betaOutUpd =
1170  betaOut +
1171  alpaka::math::copysign(
1172  acc,
1173  alpaka::math::asin(
1174  acc,
1175  alpaka::math::min(acc, sdOut_dr * k2Rinv1GeVf / alpaka::math::abs(acc, pt_betaIn), kSinAlphaMax)),
1176  betaIn); //FIXME: need a faster version
1177  betaAv = (alpaka::math::abs(acc, betaOut) > 0.2f * alpaka::math::abs(acc, betaIn))
1178  ? (0.5f * (betaInUpd + betaOutUpd))
1179  : betaInUpd;
1180 
1181  //1st update
1182  pt_beta = dr * k2Rinv1GeVf / alpaka::math::sin(acc, betaAv); //get a better pt estimate
1183  betaIn += alpaka::math::copysign(
1184  acc,
1185  alpaka::math::asin(
1186  acc, alpaka::math::min(acc, sdIn_dr * k2Rinv1GeVf / alpaka::math::abs(acc, pt_beta), kSinAlphaMax)),
1187  betaIn); //FIXME: need a faster version
1188  betaOut += alpaka::math::copysign(
1189  acc,
1190  alpaka::math::asin(
1191  acc, alpaka::math::min(acc, sdOut_dr * k2Rinv1GeVf / alpaka::math::abs(acc, pt_beta), kSinAlphaMax)),
1192  betaIn); //FIXME: need a faster version
1193  //update the av and pt
1194  betaAv = 0.5f * (betaIn + betaOut);
1195  //2nd update
1196  pt_beta = dr * k2Rinv1GeVf / alpaka::math::sin(acc, betaAv); //get a better pt estimate
1197  }
1198  }
1199 
1200  template <typename TAcc>
1201  ALPAKA_FN_ACC ALPAKA_FN_INLINE bool runQuintupletDefaultAlgoBBBB(TAcc const& acc,
1202  ModulesConst modules,
1203  MiniDoubletsConst mds,
1204  SegmentsConst segments,
1205  uint16_t innerInnerLowerModuleIndex,
1206  uint16_t innerOuterLowerModuleIndex,
1207  uint16_t outerInnerLowerModuleIndex,
1208  uint16_t outerOuterLowerModuleIndex,
1209  unsigned int innerSegmentIndex,
1210  unsigned int outerSegmentIndex,
1211  unsigned int firstMDIndex,
1212  unsigned int secondMDIndex,
1213  unsigned int thirdMDIndex,
1214  unsigned int fourthMDIndex) {
1215  bool isPS_InLo = (modules.moduleType()[innerInnerLowerModuleIndex] == PS);
1216  bool isPS_OutLo = (modules.moduleType()[outerInnerLowerModuleIndex] == PS);
1217 
1218  float rt_InLo = mds.anchorRt()[firstMDIndex];
1219  float rt_InOut = mds.anchorRt()[secondMDIndex];
1220  float rt_OutLo = mds.anchorRt()[thirdMDIndex];
1221 
1222  float z_InLo = mds.anchorZ()[firstMDIndex];
1223  float z_InOut = mds.anchorZ()[secondMDIndex];
1224  float z_OutLo = mds.anchorZ()[thirdMDIndex];
1225 
1226  float alpha1GeV_OutLo =
1227  alpaka::math::asin(acc, alpaka::math::min(acc, rt_OutLo * k2Rinv1GeVf / ptCut, kSinAlphaMax));
1228 
1229  float rtRatio_OutLoInLo = rt_OutLo / rt_InLo; // Outer segment beginning rt divided by inner segment beginning rt;
1230  float dzDrtScale =
1231  alpaka::math::tan(acc, alpha1GeV_OutLo) / alpha1GeV_OutLo; // The track can bend in r-z plane slightly
1232  float zpitch_InLo = (isPS_InLo ? kPixelPSZpitch : kStrip2SZpitch);
1233  float zpitch_OutLo = (isPS_OutLo ? kPixelPSZpitch : kStrip2SZpitch);
1234 
1235  float zHi = z_InLo + (z_InLo + kDeltaZLum) * (rtRatio_OutLoInLo - 1.f) * (z_InLo < 0.f ? 1.f : dzDrtScale) +
1236  (zpitch_InLo + zpitch_OutLo);
1237  float zLo = z_InLo + (z_InLo - kDeltaZLum) * (rtRatio_OutLoInLo - 1.f) * (z_InLo > 0.f ? 1.f : dzDrtScale) -
1238  (zpitch_InLo + zpitch_OutLo);
1239 
1240  //Cut 1 - z compatibility
1241  if ((z_OutLo < zLo) || (z_OutLo > zHi))
1242  return false;
1243 
1244  float drt_OutLo_InLo = (rt_OutLo - rt_InLo);
1245  float r3_InLo = alpaka::math::sqrt(acc, z_InLo * z_InLo + rt_InLo * rt_InLo);
1246  float drt_InSeg = rt_InOut - rt_InLo;
1247  float dz_InSeg = z_InOut - z_InLo;
1248  float dr3_InSeg = alpaka::math::sqrt(acc, rt_InOut * rt_InOut + z_InOut * z_InOut) -
1249  alpaka::math::sqrt(acc, rt_InLo * rt_InLo + z_InLo * z_InLo);
1250 
1251  float coshEta = dr3_InSeg / drt_InSeg;
1252  float dzErr = (zpitch_InLo + zpitch_OutLo) * (zpitch_InLo + zpitch_OutLo) * 2.f;
1253 
1254  float thetaMuls2 = (kMulsInGeV * kMulsInGeV) * (0.1f + 0.2f * (rt_OutLo - rt_InLo) / 50.f) * (r3_InLo / rt_InLo);
1255  float muls2 = thetaMuls2 * 9.f / (ptCut * ptCut) * 16.f;
1256  dzErr += muls2 * drt_OutLo_InLo * drt_OutLo_InLo / 3.f * coshEta * coshEta;
1257  dzErr = alpaka::math::sqrt(acc, dzErr);
1258 
1259  // Constructing upper and lower bound
1260  const float dzMean = dz_InSeg / drt_InSeg * drt_OutLo_InLo;
1261  const float zWindow =
1262  dzErr / drt_InSeg * drt_OutLo_InLo +
1263  (zpitch_InLo + zpitch_OutLo); //FIXME for ptCut lower than ~0.8 need to add curv path correction
1264  float zLoPointed = z_InLo + dzMean * (z_InLo > 0.f ? 1.f : dzDrtScale) - zWindow;
1265  float zHiPointed = z_InLo + dzMean * (z_InLo < 0.f ? 1.f : dzDrtScale) + zWindow;
1266 
1267  // Cut #2: Pointed Z (Inner segment two MD points to outer segment inner MD)
1268  if ((z_OutLo < zLoPointed) || (z_OutLo > zHiPointed))
1269  return false;
1270 
1271  float pvOffset = 0.1f / rt_OutLo;
1272  float dPhiCut = alpha1GeV_OutLo + alpaka::math::sqrt(acc, muls2 + pvOffset * pvOffset);
1273 
1274  float deltaPhiPos = phi_mpi_pi(acc, mds.anchorPhi()[fourthMDIndex] - mds.anchorPhi()[secondMDIndex]);
1275  // Cut #3: FIXME:deltaPhiPos can be tighter
1276  if (alpaka::math::abs(acc, deltaPhiPos) > dPhiCut)
1277  return false;
1278 
1279  float midPointX = 0.5f * (mds.anchorX()[firstMDIndex] + mds.anchorX()[thirdMDIndex]);
1280  float midPointY = 0.5f * (mds.anchorY()[firstMDIndex] + mds.anchorY()[thirdMDIndex]);
1281  float diffX = mds.anchorX()[thirdMDIndex] - mds.anchorX()[firstMDIndex];
1282  float diffY = mds.anchorY()[thirdMDIndex] - mds.anchorY()[firstMDIndex];
1283 
1284  float dPhi = deltaPhi(acc, midPointX, midPointY, diffX, diffY);
1285 
1286  // Cut #4: deltaPhiChange
1287  if (alpaka::math::abs(acc, dPhi) > dPhiCut)
1288  return false;
1289 
1290  // First obtaining the raw betaIn and betaOut values without any correction and just purely based on the mini-doublet hit positions
1291  float alpha_InLo = __H2F(segments.dPhiChanges()[innerSegmentIndex]);
1292  float alpha_OutLo = __H2F(segments.dPhiChanges()[outerSegmentIndex]);
1293 
1294  bool isEC_lastLayer = modules.subdets()[outerOuterLowerModuleIndex] == Endcap and
1295  modules.moduleType()[outerOuterLowerModuleIndex] == TwoS;
1296 
1297  float alpha_OutUp, alpha_OutUp_highEdge, alpha_OutUp_lowEdge;
1298 
1299  alpha_OutUp = phi_mpi_pi(acc,
1300  phi(acc,
1301  mds.anchorX()[fourthMDIndex] - mds.anchorX()[thirdMDIndex],
1302  mds.anchorY()[fourthMDIndex] - mds.anchorY()[thirdMDIndex]) -
1303  mds.anchorPhi()[fourthMDIndex]);
1304 
1305  alpha_OutUp_highEdge = alpha_OutUp;
1306  alpha_OutUp_lowEdge = alpha_OutUp;
1307 
1308  float tl_axis_x = mds.anchorX()[fourthMDIndex] - mds.anchorX()[firstMDIndex];
1309  float tl_axis_y = mds.anchorY()[fourthMDIndex] - mds.anchorY()[firstMDIndex];
1310  float tl_axis_highEdge_x = tl_axis_x;
1311  float tl_axis_highEdge_y = tl_axis_y;
1312  float tl_axis_lowEdge_x = tl_axis_x;
1313  float tl_axis_lowEdge_y = tl_axis_y;
1314 
1315  float betaIn = alpha_InLo - phi_mpi_pi(acc, phi(acc, tl_axis_x, tl_axis_y) - mds.anchorPhi()[firstMDIndex]);
1316 
1317  float betaInRHmin = betaIn;
1318  float betaInRHmax = betaIn;
1319  float betaOut = -alpha_OutUp + phi_mpi_pi(acc, phi(acc, tl_axis_x, tl_axis_y) - mds.anchorPhi()[fourthMDIndex]);
1320 
1321  float betaOutRHmin = betaOut;
1322  float betaOutRHmax = betaOut;
1323 
1324  if (isEC_lastLayer) {
1325  alpha_OutUp_highEdge = phi_mpi_pi(acc,
1326  phi(acc,
1327  mds.anchorHighEdgeX()[fourthMDIndex] - mds.anchorX()[thirdMDIndex],
1328  mds.anchorHighEdgeY()[fourthMDIndex] - mds.anchorY()[thirdMDIndex]) -
1329  mds.anchorHighEdgePhi()[fourthMDIndex]);
1330  alpha_OutUp_lowEdge = phi_mpi_pi(acc,
1331  phi(acc,
1332  mds.anchorLowEdgeX()[fourthMDIndex] - mds.anchorX()[thirdMDIndex],
1333  mds.anchorLowEdgeY()[fourthMDIndex] - mds.anchorY()[thirdMDIndex]) -
1334  mds.anchorLowEdgePhi()[fourthMDIndex]);
1335 
1336  tl_axis_highEdge_x = mds.anchorHighEdgeX()[fourthMDIndex] - mds.anchorX()[firstMDIndex];
1337  tl_axis_highEdge_y = mds.anchorHighEdgeY()[fourthMDIndex] - mds.anchorY()[firstMDIndex];
1338  tl_axis_lowEdge_x = mds.anchorLowEdgeX()[fourthMDIndex] - mds.anchorX()[firstMDIndex];
1339  tl_axis_lowEdge_y = mds.anchorLowEdgeY()[fourthMDIndex] - mds.anchorY()[firstMDIndex];
1340 
1341  betaOutRHmin =
1342  -alpha_OutUp_highEdge +
1343  phi_mpi_pi(acc, phi(acc, tl_axis_highEdge_x, tl_axis_highEdge_y) - mds.anchorHighEdgePhi()[fourthMDIndex]);
1344  betaOutRHmax =
1345  -alpha_OutUp_lowEdge +
1346  phi_mpi_pi(acc, phi(acc, tl_axis_lowEdge_x, tl_axis_lowEdge_y) - mds.anchorLowEdgePhi()[fourthMDIndex]);
1347  }
1348 
1349  //beta computation
1350  float drt_tl_axis = alpaka::math::sqrt(acc, tl_axis_x * tl_axis_x + tl_axis_y * tl_axis_y);
1351 
1352  float corrF = 1.f;
1353  //innerOuterAnchor - innerInnerAnchor
1354  const float rt_InSeg = alpaka::math::sqrt(acc,
1355  (mds.anchorX()[secondMDIndex] - mds.anchorX()[firstMDIndex]) *
1356  (mds.anchorX()[secondMDIndex] - mds.anchorX()[firstMDIndex]) +
1357  (mds.anchorY()[secondMDIndex] - mds.anchorY()[firstMDIndex]) *
1358  (mds.anchorY()[secondMDIndex] - mds.anchorY()[firstMDIndex]));
1359  float betaInCut =
1360  alpaka::math::asin(
1361  acc, alpaka::math::min(acc, (-rt_InSeg * corrF + drt_tl_axis) * k2Rinv1GeVf / ptCut, kSinAlphaMax)) +
1362  (0.02f / drt_InSeg);
1363 
1364  //Cut #5: first beta cut
1365  if (alpaka::math::abs(acc, betaInRHmin) >= betaInCut)
1366  return false;
1367 
1368  float betaAv = 0.5f * (betaIn + betaOut);
1369  float pt_beta = drt_tl_axis * k2Rinv1GeVf / alpaka::math::sin(acc, betaAv);
1370  int lIn = 5;
1371  int lOut = isEC_lastLayer ? 11 : 5;
1372  float sdOut_dr = alpaka::math::sqrt(acc,
1373  (mds.anchorX()[fourthMDIndex] - mds.anchorX()[thirdMDIndex]) *
1374  (mds.anchorX()[fourthMDIndex] - mds.anchorX()[thirdMDIndex]) +
1375  (mds.anchorY()[fourthMDIndex] - mds.anchorY()[thirdMDIndex]) *
1376  (mds.anchorY()[fourthMDIndex] - mds.anchorY()[thirdMDIndex]));
1377  float sdOut_d = mds.anchorRt()[fourthMDIndex] - mds.anchorRt()[thirdMDIndex];
1378 
1379  runDeltaBetaIterationsT5(acc, betaIn, betaOut, betaAv, pt_beta, rt_InSeg, sdOut_dr, drt_tl_axis, lIn);
1380 
1381  const float betaInMMSF = (alpaka::math::abs(acc, betaInRHmin + betaInRHmax) > 0)
1382  ? (2.f * betaIn / alpaka::math::abs(acc, betaInRHmin + betaInRHmax))
1383  : 0.f; //mean value of min,max is the old betaIn
1384  const float betaOutMMSF = (alpaka::math::abs(acc, betaOutRHmin + betaOutRHmax) > 0)
1385  ? (2.f * betaOut / alpaka::math::abs(acc, betaOutRHmin + betaOutRHmax))
1386  : 0.f;
1387  betaInRHmin *= betaInMMSF;
1388  betaInRHmax *= betaInMMSF;
1389  betaOutRHmin *= betaOutMMSF;
1390  betaOutRHmax *= betaOutMMSF;
1391 
1392  float min_ptBeta_maxPtBeta = alpaka::math::min(
1393  acc, alpaka::math::abs(acc, pt_beta), kPt_betaMax); //need to confimm the range-out value of 7 GeV
1394  const float dBetaMuls2 = thetaMuls2 * 16.f / (min_ptBeta_maxPtBeta * min_ptBeta_maxPtBeta);
1395 
1396  const float alphaInAbsReg =
1397  alpaka::math::max(acc,
1398  alpaka::math::abs(acc, alpha_InLo),
1399  alpaka::math::asin(acc, alpaka::math::min(acc, rt_InLo * k2Rinv1GeVf / 3.0f, kSinAlphaMax)));
1400  const float alphaOutAbsReg =
1401  alpaka::math::max(acc,
1402  alpaka::math::abs(acc, alpha_OutLo),
1403  alpaka::math::asin(acc, alpaka::math::min(acc, rt_OutLo * k2Rinv1GeVf / 3.0f, kSinAlphaMax)));
1404  const float dBetaInLum = lIn < 11 ? 0.0f : alpaka::math::abs(acc, alphaInAbsReg * kDeltaZLum / z_InLo);
1405  const float dBetaOutLum = lOut < 11 ? 0.0f : alpaka::math::abs(acc, alphaOutAbsReg * kDeltaZLum / z_OutLo);
1406  const float dBetaLum2 = (dBetaInLum + dBetaOutLum) * (dBetaInLum + dBetaOutLum);
1407  const float sinDPhi = alpaka::math::sin(acc, dPhi);
1408 
1409  const float dBetaRIn2 = 0; // TODO-RH
1410  float dBetaROut = 0;
1411  if (isEC_lastLayer) {
1412  dBetaROut = (alpaka::math::sqrt(acc,
1413  mds.anchorHighEdgeX()[fourthMDIndex] * mds.anchorHighEdgeX()[fourthMDIndex] +
1414  mds.anchorHighEdgeY()[fourthMDIndex] * mds.anchorHighEdgeY()[fourthMDIndex]) -
1415  alpaka::math::sqrt(acc,
1416  mds.anchorLowEdgeX()[fourthMDIndex] * mds.anchorLowEdgeX()[fourthMDIndex] +
1417  mds.anchorLowEdgeY()[fourthMDIndex] * mds.anchorLowEdgeY()[fourthMDIndex])) *
1418  sinDPhi / drt_tl_axis;
1419  }
1420 
1421  const float dBetaROut2 = dBetaROut * dBetaROut;
1422 
1423  float betaOutCut =
1424  alpaka::math::asin(acc, alpaka::math::min(acc, drt_tl_axis * k2Rinv1GeVf / ptCut, kSinAlphaMax)) +
1425  (0.02f / sdOut_d) + alpaka::math::sqrt(acc, dBetaLum2 + dBetaMuls2);
1426 
1427  //Cut #6: The real beta cut
1428  if (alpaka::math::abs(acc, betaOut) >= betaOutCut)
1429  return false;
1430 
1431  float dBetaRes = 0.02f / alpaka::math::min(acc, sdOut_d, drt_InSeg);
1432  float dBetaCut2 =
1433  (dBetaRes * dBetaRes * 2.0f + dBetaMuls2 + dBetaLum2 + dBetaRIn2 + dBetaROut2 +
1434  0.25f *
1435  (alpaka::math::abs(acc, betaInRHmin - betaInRHmax) + alpaka::math::abs(acc, betaOutRHmin - betaOutRHmax)) *
1436  (alpaka::math::abs(acc, betaInRHmin - betaInRHmax) + alpaka::math::abs(acc, betaOutRHmin - betaOutRHmax)));
1437 
1438  float dBeta = betaIn - betaOut;
1439  return dBeta * dBeta <= dBetaCut2;
1440  }
1441 
1442  template <typename TAcc>
1443  ALPAKA_FN_ACC ALPAKA_FN_INLINE bool runQuintupletDefaultAlgoBBEE(TAcc const& acc,
1444  ModulesConst modules,
1445  MiniDoubletsConst mds,
1446  SegmentsConst segments,
1447  uint16_t innerInnerLowerModuleIndex,
1448  uint16_t innerOuterLowerModuleIndex,
1449  uint16_t outerInnerLowerModuleIndex,
1450  uint16_t outerOuterLowerModuleIndex,
1451  unsigned int innerSegmentIndex,
1452  unsigned int outerSegmentIndex,
1453  unsigned int firstMDIndex,
1454  unsigned int secondMDIndex,
1455  unsigned int thirdMDIndex,
1456  unsigned int fourthMDIndex) {
1457  bool isPS_InLo = (modules.moduleType()[innerInnerLowerModuleIndex] == PS);
1458  bool isPS_OutLo = (modules.moduleType()[outerInnerLowerModuleIndex] == PS);
1459 
1460  float rt_InLo = mds.anchorRt()[firstMDIndex];
1461  float rt_InOut = mds.anchorRt()[secondMDIndex];
1462  float rt_OutLo = mds.anchorRt()[thirdMDIndex];
1463 
1464  float z_InLo = mds.anchorZ()[firstMDIndex];
1465  float z_InOut = mds.anchorZ()[secondMDIndex];
1466  float z_OutLo = mds.anchorZ()[thirdMDIndex];
1467 
1468  float alpha1GeV_OutLo =
1469  alpaka::math::asin(acc, alpaka::math::min(acc, rt_OutLo * k2Rinv1GeVf / ptCut, kSinAlphaMax));
1470 
1471  float dzDrtScale =
1472  alpaka::math::tan(acc, alpha1GeV_OutLo) / alpha1GeV_OutLo; // The track can bend in r-z plane slightly
1473  float zpitch_InLo = (isPS_InLo ? kPixelPSZpitch : kStrip2SZpitch);
1474  float zpitch_OutLo = (isPS_OutLo ? kPixelPSZpitch : kStrip2SZpitch);
1475  float zGeom = zpitch_InLo + zpitch_OutLo;
1476 
1477  // Cut #0: Preliminary (Only here in endcap case)
1478  if (z_InLo * z_OutLo <= 0)
1479  return false;
1480 
1481  float dLum = alpaka::math::copysign(acc, kDeltaZLum, z_InLo);
1482  bool isOutSgInnerMDPS = modules.moduleType()[outerInnerLowerModuleIndex] == PS;
1483  float rtGeom1 = isOutSgInnerMDPS ? kPixelPSZpitch : kStrip2SZpitch;
1484  float zGeom1 = alpaka::math::copysign(acc, zGeom, z_InLo);
1485  float rtLo = rt_InLo * (1.f + (z_OutLo - z_InLo - zGeom1) / (z_InLo + zGeom1 + dLum) / dzDrtScale) -
1486  rtGeom1; //slope correction only on the lower end
1487  float rtOut = rt_OutLo;
1488 
1489  //Cut #1: rt condition
1490  if (rtOut < rtLo)
1491  return false;
1492 
1493  float zInForHi = z_InLo - zGeom1 - dLum;
1494  if (zInForHi * z_InLo < 0) {
1495  zInForHi = alpaka::math::copysign(acc, 0.1f, z_InLo);
1496  }
1497  float rtHi = rt_InLo * (1.f + (z_OutLo - z_InLo + zGeom1) / zInForHi) + rtGeom1;
1498 
1499  //Cut #2: rt condition
1500  if ((rt_OutLo < rtLo) || (rt_OutLo > rtHi))
1501  return false;
1502 
1503  float rIn = alpaka::math::sqrt(acc, z_InLo * z_InLo + rt_InLo * rt_InLo);
1504  const float drtSDIn = rt_InOut - rt_InLo;
1505  const float dzSDIn = z_InOut - z_InLo;
1506  const float dr3SDIn = alpaka::math::sqrt(acc, rt_InOut * rt_InOut + z_InOut * z_InOut) -
1507  alpaka::math::sqrt(acc, rt_InLo * rt_InLo + z_InLo * z_InLo);
1508 
1509  const float coshEta = dr3SDIn / drtSDIn; //direction estimate
1510  const float dzOutInAbs = alpaka::math::abs(acc, z_OutLo - z_InLo);
1511  const float multDzDr = dzOutInAbs * coshEta / (coshEta * coshEta - 1.f);
1512  const float zGeom1_another = kPixelPSZpitch;
1513  float kZ = (z_OutLo - z_InLo) / dzSDIn;
1514  float drtErr =
1515  zGeom1_another * zGeom1_another * drtSDIn * drtSDIn / dzSDIn / dzSDIn * (1.f - 2.f * kZ + 2.f * kZ * kZ);
1516  const float thetaMuls2 = (kMulsInGeV * kMulsInGeV) * (0.1f + 0.2f * (rt_OutLo - rt_InLo) / 50.f) * (rIn / rt_InLo);
1517  const float muls2 = thetaMuls2 * 9.f / (ptCut * ptCut) * 16.f;
1518  drtErr += muls2 * multDzDr * multDzDr / 3.f * coshEta * coshEta;
1519  drtErr = alpaka::math::sqrt(acc, drtErr);
1520 
1521  //Cut #3: rt-z pointed
1522  if ((kZ < 0) || (rtOut < rtLo) || (rtOut > rtHi))
1523  return false;
1524 
1525  const float pvOffset = 0.1f / rt_OutLo;
1526  float dPhiCut = alpha1GeV_OutLo + alpaka::math::sqrt(acc, muls2 + pvOffset * pvOffset);
1527 
1528  float deltaPhiPos = phi_mpi_pi(acc, mds.anchorPhi()[fourthMDIndex] - mds.anchorPhi()[secondMDIndex]);
1529 
1530  //Cut #4: deltaPhiPos can be tighter
1531  if (alpaka::math::abs(acc, deltaPhiPos) > dPhiCut)
1532  return false;
1533 
1534  float midPointX = 0.5f * (mds.anchorX()[firstMDIndex] + mds.anchorX()[thirdMDIndex]);
1535  float midPointY = 0.5f * (mds.anchorY()[firstMDIndex] + mds.anchorY()[thirdMDIndex]);
1536  float diffX = mds.anchorX()[thirdMDIndex] - mds.anchorX()[firstMDIndex];
1537  float diffY = mds.anchorY()[thirdMDIndex] - mds.anchorY()[firstMDIndex];
1538 
1539  float dPhi = deltaPhi(acc, midPointX, midPointY, diffX, diffY);
1540  // Cut #5: deltaPhiChange
1541  if (alpaka::math::abs(acc, dPhi) > dPhiCut)
1542  return false;
1543 
1544  float sdIn_alpha = __H2F(segments.dPhiChanges()[innerSegmentIndex]);
1545  float sdIn_alpha_min = __H2F(segments.dPhiChangeMins()[innerSegmentIndex]);
1546  float sdIn_alpha_max = __H2F(segments.dPhiChangeMaxs()[innerSegmentIndex]);
1547  float sdOut_alpha = sdIn_alpha;
1548 
1549  float sdOut_alphaOut = phi_mpi_pi(acc,
1550  phi(acc,
1551  mds.anchorX()[fourthMDIndex] - mds.anchorX()[thirdMDIndex],
1552  mds.anchorY()[fourthMDIndex] - mds.anchorY()[thirdMDIndex]) -
1553  mds.anchorPhi()[fourthMDIndex]);
1554 
1555  float sdOut_alphaOut_min = phi_mpi_pi(
1556  acc, __H2F(segments.dPhiChangeMins()[outerSegmentIndex]) - __H2F(segments.dPhiMins()[outerSegmentIndex]));
1557  float sdOut_alphaOut_max = phi_mpi_pi(
1558  acc, __H2F(segments.dPhiChangeMaxs()[outerSegmentIndex]) - __H2F(segments.dPhiMaxs()[outerSegmentIndex]));
1559 
1560  float tl_axis_x = mds.anchorX()[fourthMDIndex] - mds.anchorX()[firstMDIndex];
1561  float tl_axis_y = mds.anchorY()[fourthMDIndex] - mds.anchorY()[firstMDIndex];
1562 
1563  float betaIn = sdIn_alpha - phi_mpi_pi(acc, phi(acc, tl_axis_x, tl_axis_y) - mds.anchorPhi()[firstMDIndex]);
1564 
1565  float betaInRHmin = betaIn;
1566  float betaInRHmax = betaIn;
1567  float betaOut = -sdOut_alphaOut + phi_mpi_pi(acc, phi(acc, tl_axis_x, tl_axis_y) - mds.anchorPhi()[fourthMDIndex]);
1568 
1569  float betaOutRHmin = betaOut;
1570  float betaOutRHmax = betaOut;
1571 
1572  bool isEC_secondLayer = (modules.subdets()[innerOuterLowerModuleIndex] == Endcap) and
1573  (modules.moduleType()[innerOuterLowerModuleIndex] == TwoS);
1574 
1575  if (isEC_secondLayer) {
1576  betaInRHmin = betaIn - sdIn_alpha_min + sdIn_alpha;
1577  betaInRHmax = betaIn - sdIn_alpha_max + sdIn_alpha;
1578  }
1579 
1580  betaOutRHmin = betaOut - sdOut_alphaOut_min + sdOut_alphaOut;
1581  betaOutRHmax = betaOut - sdOut_alphaOut_max + sdOut_alphaOut;
1582 
1583  float swapTemp;
1584  if (alpaka::math::abs(acc, betaOutRHmin) > alpaka::math::abs(acc, betaOutRHmax)) {
1585  swapTemp = betaOutRHmin;
1586  betaOutRHmin = betaOutRHmax;
1587  betaOutRHmax = swapTemp;
1588  }
1589 
1590  if (alpaka::math::abs(acc, betaInRHmin) > alpaka::math::abs(acc, betaInRHmax)) {
1591  swapTemp = betaInRHmin;
1592  betaInRHmin = betaInRHmax;
1593  betaInRHmax = swapTemp;
1594  }
1595 
1596  float sdIn_dr = alpaka::math::sqrt(acc,
1597  (mds.anchorX()[secondMDIndex] - mds.anchorX()[firstMDIndex]) *
1598  (mds.anchorX()[secondMDIndex] - mds.anchorX()[firstMDIndex]) +
1599  (mds.anchorY()[secondMDIndex] - mds.anchorY()[firstMDIndex]) *
1600  (mds.anchorY()[secondMDIndex] - mds.anchorY()[firstMDIndex]));
1601  float sdIn_d = rt_InOut - rt_InLo;
1602 
1603  float dr = alpaka::math::sqrt(acc, tl_axis_x * tl_axis_x + tl_axis_y * tl_axis_y);
1604  const float corrF = 1.f;
1605  float betaInCut =
1606  alpaka::math::asin(acc, alpaka::math::min(acc, (-sdIn_dr * corrF + dr) * k2Rinv1GeVf / ptCut, kSinAlphaMax)) +
1607  (0.02f / sdIn_d);
1608 
1609  //Cut #6: first beta cut
1610  if (alpaka::math::abs(acc, betaInRHmin) >= betaInCut)
1611  return false;
1612 
1613  float betaAv = 0.5f * (betaIn + betaOut);
1614  float pt_beta = dr * k2Rinv1GeVf / alpaka::math::sin(acc, betaAv);
1615 
1616  float lIn = 5;
1617  float lOut = 11;
1618 
1619  float sdOut_dr = alpaka::math::sqrt(acc,
1620  (mds.anchorX()[fourthMDIndex] - mds.anchorX()[thirdMDIndex]) *
1621  (mds.anchorX()[fourthMDIndex] - mds.anchorX()[thirdMDIndex]) +
1622  (mds.anchorY()[fourthMDIndex] - mds.anchorY()[thirdMDIndex]) *
1623  (mds.anchorY()[fourthMDIndex] - mds.anchorY()[thirdMDIndex]));
1624  float sdOut_d = mds.anchorRt()[fourthMDIndex] - mds.anchorRt()[thirdMDIndex];
1625 
1626  runDeltaBetaIterationsT5(acc, betaIn, betaOut, betaAv, pt_beta, sdIn_dr, sdOut_dr, dr, lIn);
1627 
1628  const float betaInMMSF = (alpaka::math::abs(acc, betaInRHmin + betaInRHmax) > 0)
1629  ? (2.f * betaIn / alpaka::math::abs(acc, betaInRHmin + betaInRHmax))
1630  : 0.; //mean value of min,max is the old betaIn
1631  const float betaOutMMSF = (alpaka::math::abs(acc, betaOutRHmin + betaOutRHmax) > 0)
1632  ? (2.f * betaOut / alpaka::math::abs(acc, betaOutRHmin + betaOutRHmax))
1633  : 0.;
1634  betaInRHmin *= betaInMMSF;
1635  betaInRHmax *= betaInMMSF;
1636  betaOutRHmin *= betaOutMMSF;
1637  betaOutRHmax *= betaOutMMSF;
1638 
1639  float min_ptBeta_maxPtBeta = alpaka::math::min(
1640  acc, alpaka::math::abs(acc, pt_beta), kPt_betaMax); //need to confirm the range-out value of 7 GeV
1641  const float dBetaMuls2 = thetaMuls2 * 16.f / (min_ptBeta_maxPtBeta * min_ptBeta_maxPtBeta);
1642 
1643  const float alphaInAbsReg =
1644  alpaka::math::max(acc,
1645  alpaka::math::abs(acc, sdIn_alpha),
1646  alpaka::math::asin(acc, alpaka::math::min(acc, rt_InLo * k2Rinv1GeVf / 3.0f, kSinAlphaMax)));
1647  const float alphaOutAbsReg =
1648  alpaka::math::max(acc,
1649  alpaka::math::abs(acc, sdOut_alpha),
1650  alpaka::math::asin(acc, alpaka::math::min(acc, rt_OutLo * k2Rinv1GeVf / 3.0f, kSinAlphaMax)));
1651  const float dBetaInLum = lIn < 11 ? 0.0f : alpaka::math::abs(acc, alphaInAbsReg * kDeltaZLum / z_InLo);
1652  const float dBetaOutLum = lOut < 11 ? 0.0f : alpaka::math::abs(acc, alphaOutAbsReg * kDeltaZLum / z_OutLo);
1653  const float dBetaLum2 = (dBetaInLum + dBetaOutLum) * (dBetaInLum + dBetaOutLum);
1654  const float sinDPhi = alpaka::math::sin(acc, dPhi);
1655 
1656  const float dBetaRIn2 = 0; // TODO-RH
1657  float dBetaROut = 0;
1658  if (modules.moduleType()[outerOuterLowerModuleIndex] == TwoS) {
1659  dBetaROut = (alpaka::math::sqrt(acc,
1660  mds.anchorHighEdgeX()[fourthMDIndex] * mds.anchorHighEdgeX()[fourthMDIndex] +
1661  mds.anchorHighEdgeY()[fourthMDIndex] * mds.anchorHighEdgeY()[fourthMDIndex]) -
1662  alpaka::math::sqrt(acc,
1663  mds.anchorLowEdgeX()[fourthMDIndex] * mds.anchorLowEdgeX()[fourthMDIndex] +
1664  mds.anchorLowEdgeY()[fourthMDIndex] * mds.anchorLowEdgeY()[fourthMDIndex])) *
1665  sinDPhi / dr;
1666  }
1667 
1668  const float dBetaROut2 = dBetaROut * dBetaROut;
1669  float betaOutCut = alpaka::math::asin(acc, alpaka::math::min(acc, dr * k2Rinv1GeVf / ptCut, kSinAlphaMax)) +
1670  (0.02f / sdOut_d) + alpaka::math::sqrt(acc, dBetaLum2 + dBetaMuls2);
1671 
1672  //Cut #6: The real beta cut
1673  if (alpaka::math::abs(acc, betaOut) >= betaOutCut)
1674  return false;
1675 
1676  float dBetaRes = 0.02f / alpaka::math::min(acc, sdOut_d, sdIn_d);
1677  float dBetaCut2 =
1678  (dBetaRes * dBetaRes * 2.0f + dBetaMuls2 + dBetaLum2 + dBetaRIn2 + dBetaROut2 +
1679  0.25f *
1680  (alpaka::math::abs(acc, betaInRHmin - betaInRHmax) + alpaka::math::abs(acc, betaOutRHmin - betaOutRHmax)) *
1681  (alpaka::math::abs(acc, betaInRHmin - betaInRHmax) + alpaka::math::abs(acc, betaOutRHmin - betaOutRHmax)));
1682  float dBeta = betaIn - betaOut;
1683  //Cut #7: Cut on dBet
1684  return dBeta * dBeta <= dBetaCut2;
1685  }
1686 
1687  template <typename TAcc>
1688  ALPAKA_FN_ACC ALPAKA_FN_INLINE bool runQuintupletDefaultAlgoEEEE(TAcc const& acc,
1689  ModulesConst modules,
1690  MiniDoubletsConst mds,
1691  SegmentsConst segments,
1692  uint16_t innerInnerLowerModuleIndex,
1693  uint16_t innerOuterLowerModuleIndex,
1694  uint16_t outerInnerLowerModuleIndex,
1695  uint16_t outerOuterLowerModuleIndex,
1696  unsigned int innerSegmentIndex,
1697  unsigned int outerSegmentIndex,
1698  unsigned int firstMDIndex,
1699  unsigned int secondMDIndex,
1700  unsigned int thirdMDIndex,
1701  unsigned int fourthMDIndex) {
1702  float rt_InLo = mds.anchorRt()[firstMDIndex];
1703  float rt_InOut = mds.anchorRt()[secondMDIndex];
1704  float rt_OutLo = mds.anchorRt()[thirdMDIndex];
1705 
1706  float z_InLo = mds.anchorZ()[firstMDIndex];
1707  float z_InOut = mds.anchorZ()[secondMDIndex];
1708  float z_OutLo = mds.anchorZ()[thirdMDIndex];
1709 
1710  float alpha1GeV_OutLo =
1711  alpaka::math::asin(acc, alpaka::math::min(acc, rt_OutLo * k2Rinv1GeVf / ptCut, kSinAlphaMax));
1712 
1713  float dzDrtScale =
1714  alpaka::math::tan(acc, alpha1GeV_OutLo) / alpha1GeV_OutLo; // The track can bend in r-z plane slightly
1715 
1716  // Cut #0: Preliminary (Only here in endcap case)
1717  if ((z_InLo * z_OutLo) <= 0)
1718  return false;
1719 
1720  float dLum = alpaka::math::copysign(acc, kDeltaZLum, z_InLo);
1721  bool isOutSgInnerMDPS = modules.moduleType()[outerInnerLowerModuleIndex] == PS;
1722  bool isInSgInnerMDPS = modules.moduleType()[innerInnerLowerModuleIndex] == PS;
1723 
1724  float rtGeom = (isInSgInnerMDPS and isOutSgInnerMDPS) ? 2.f * kPixelPSZpitch
1725  : (isInSgInnerMDPS or isOutSgInnerMDPS) ? kPixelPSZpitch + kStrip2SZpitch
1726  : 2.f * kStrip2SZpitch;
1727 
1728  float dz = z_OutLo - z_InLo;
1729  float rtLo = rt_InLo * (1.f + dz / (z_InLo + dLum) / dzDrtScale) - rtGeom; //slope correction only on the lower end
1730 
1731  float rtOut = rt_OutLo;
1732 
1733  //Cut #1: rt condition
1734 
1735  float rtHi = rt_InLo * (1.f + dz / (z_InLo - dLum)) + rtGeom;
1736 
1737  if ((rtOut < rtLo) || (rtOut > rtHi))
1738  return false;
1739 
1740  bool isInSgOuterMDPS = modules.moduleType()[innerOuterLowerModuleIndex] == PS;
1741 
1742  const float drtSDIn = rt_InOut - rt_InLo;
1743  const float dzSDIn = z_InOut - z_InLo;
1744  const float dr3SDIn = alpaka::math::sqrt(acc, rt_InOut * rt_InOut + z_InOut * z_InOut) -
1745  alpaka::math::sqrt(acc, rt_InLo * rt_InLo + z_InLo * z_InLo);
1746  float coshEta = dr3SDIn / drtSDIn; //direction estimate
1747  float dzOutInAbs = alpaka::math::abs(acc, z_OutLo - z_InLo);
1748  float multDzDr = dzOutInAbs * coshEta / (coshEta * coshEta - 1.f);
1749 
1750  float kZ = (z_OutLo - z_InLo) / dzSDIn;
1751  float thetaMuls2 = (kMulsInGeV * kMulsInGeV) * (0.1f + 0.2f * (rt_OutLo - rt_InLo) / 50.f);
1752 
1753  float muls2 = thetaMuls2 * 9.f / (ptCut * ptCut) * 16.f;
1754 
1755  float drtErr =
1756  alpaka::math::sqrt(acc,
1757  kPixelPSZpitch * kPixelPSZpitch * 2.f / (dzSDIn * dzSDIn) * (dzOutInAbs * dzOutInAbs) +
1758  muls2 * multDzDr * multDzDr / 3.f * coshEta * coshEta);
1759 
1760  float drtMean = drtSDIn * dzOutInAbs / alpaka::math::abs(acc, dzSDIn);
1761  float rtWindow = drtErr + rtGeom;
1762  float rtLo_point = rt_InLo + drtMean / dzDrtScale - rtWindow;
1763  float rtHi_point = rt_InLo + drtMean + rtWindow;
1764 
1765  // Cut #3: rt-z pointed
1766  // https://github.com/slava77/cms-tkph2-ntuple/blob/superDoubletLinked-91X-noMock/doubletAnalysis.C#L3765
1767 
1768  if (isInSgInnerMDPS and isInSgOuterMDPS) // If both PS then we can point
1769  {
1770  if (kZ < 0 || rtOut < rtLo_point || rtOut > rtHi_point)
1771  return false;
1772  }
1773 
1774  float pvOffset = 0.1f / rtOut;
1775  float dPhiCut = alpha1GeV_OutLo + alpaka::math::sqrt(acc, muls2 + pvOffset * pvOffset);
1776 
1777  float deltaPhiPos = phi_mpi_pi(acc, mds.anchorPhi()[fourthMDIndex] - mds.anchorPhi()[secondMDIndex]);
1778 
1779  if (alpaka::math::abs(acc, deltaPhiPos) > dPhiCut)
1780  return false;
1781 
1782  float midPointX = 0.5f * (mds.anchorX()[firstMDIndex] + mds.anchorX()[thirdMDIndex]);
1783  float midPointY = 0.5f * (mds.anchorY()[firstMDIndex] + mds.anchorY()[thirdMDIndex]);
1784  float diffX = mds.anchorX()[thirdMDIndex] - mds.anchorX()[firstMDIndex];
1785  float diffY = mds.anchorY()[thirdMDIndex] - mds.anchorY()[firstMDIndex];
1786 
1787  float dPhi = deltaPhi(acc, midPointX, midPointY, diffX, diffY);
1788 
1789  // Cut #5: deltaPhiChange
1790  if (alpaka::math::abs(acc, dPhi) > dPhiCut)
1791  return false;
1792 
1793  float sdIn_alpha = __H2F(segments.dPhiChanges()[innerSegmentIndex]);
1794  float sdOut_alpha = sdIn_alpha; //weird
1795  float sdOut_dPhiPos = phi_mpi_pi(acc, mds.anchorPhi()[fourthMDIndex] - mds.anchorPhi()[thirdMDIndex]);
1796 
1797  float sdOut_dPhiChange = __H2F(segments.dPhiChanges()[outerSegmentIndex]);
1798  float sdOut_dPhiChange_min = __H2F(segments.dPhiChangeMins()[outerSegmentIndex]);
1799  float sdOut_dPhiChange_max = __H2F(segments.dPhiChangeMaxs()[outerSegmentIndex]);
1800 
1801  float sdOut_alphaOutRHmin = phi_mpi_pi(acc, sdOut_dPhiChange_min - sdOut_dPhiPos);
1802  float sdOut_alphaOutRHmax = phi_mpi_pi(acc, sdOut_dPhiChange_max - sdOut_dPhiPos);
1803  float sdOut_alphaOut = phi_mpi_pi(acc, sdOut_dPhiChange - sdOut_dPhiPos);
1804 
1805  float tl_axis_x = mds.anchorX()[fourthMDIndex] - mds.anchorX()[firstMDIndex];
1806  float tl_axis_y = mds.anchorY()[fourthMDIndex] - mds.anchorY()[firstMDIndex];
1807 
1808  float betaIn = sdIn_alpha - phi_mpi_pi(acc, phi(acc, tl_axis_x, tl_axis_y) - mds.anchorPhi()[firstMDIndex]);
1809 
1810  float sdIn_alphaRHmin = __H2F(segments.dPhiChangeMins()[innerSegmentIndex]);
1811  float sdIn_alphaRHmax = __H2F(segments.dPhiChangeMaxs()[innerSegmentIndex]);
1812  float betaInRHmin = betaIn + sdIn_alphaRHmin - sdIn_alpha;
1813  float betaInRHmax = betaIn + sdIn_alphaRHmax - sdIn_alpha;
1814 
1815  float betaOut = -sdOut_alphaOut + phi_mpi_pi(acc, phi(acc, tl_axis_x, tl_axis_y) - mds.anchorPhi()[fourthMDIndex]);
1816 
1817  float betaOutRHmin = betaOut - sdOut_alphaOutRHmin + sdOut_alphaOut;
1818  float betaOutRHmax = betaOut - sdOut_alphaOutRHmax + sdOut_alphaOut;
1819 
1820  float swapTemp;
1821  if (alpaka::math::abs(acc, betaOutRHmin) > alpaka::math::abs(acc, betaOutRHmax)) {
1822  swapTemp = betaOutRHmin;
1823  betaOutRHmin = betaOutRHmax;
1824  betaOutRHmax = swapTemp;
1825  }
1826 
1827  if (alpaka::math::abs(acc, betaInRHmin) > alpaka::math::abs(acc, betaInRHmax)) {
1828  swapTemp = betaInRHmin;
1829  betaInRHmin = betaInRHmax;
1830  betaInRHmax = swapTemp;
1831  }
1832  float sdIn_dr = alpaka::math::sqrt(acc,
1833  (mds.anchorX()[secondMDIndex] - mds.anchorX()[firstMDIndex]) *
1834  (mds.anchorX()[secondMDIndex] - mds.anchorX()[firstMDIndex]) +
1835  (mds.anchorY()[secondMDIndex] - mds.anchorY()[firstMDIndex]) *
1836  (mds.anchorY()[secondMDIndex] - mds.anchorY()[firstMDIndex]));
1837  float sdIn_d = rt_InOut - rt_InLo;
1838 
1839  float dr = alpaka::math::sqrt(acc, tl_axis_x * tl_axis_x + tl_axis_y * tl_axis_y);
1840  const float corrF = 1.f;
1841  float betaInCut =
1842  alpaka::math::asin(acc, alpaka::math::min(acc, (-sdIn_dr * corrF + dr) * k2Rinv1GeVf / ptCut, kSinAlphaMax)) +
1843  (0.02f / sdIn_d);
1844 
1845  //Cut #6: first beta cut
1846  if (alpaka::math::abs(acc, betaInRHmin) >= betaInCut)
1847  return false;
1848 
1849  float betaAv = 0.5f * (betaIn + betaOut);
1850  float pt_beta = dr * k2Rinv1GeVf / alpaka::math::sin(acc, betaAv);
1851 
1852  int lIn = 11; //endcap
1853  int lOut = 13; //endcap
1854 
1855  float sdOut_dr = alpaka::math::sqrt(acc,
1856  (mds.anchorX()[fourthMDIndex] - mds.anchorX()[thirdMDIndex]) *
1857  (mds.anchorX()[fourthMDIndex] - mds.anchorX()[thirdMDIndex]) +
1858  (mds.anchorY()[fourthMDIndex] - mds.anchorY()[thirdMDIndex]) *
1859  (mds.anchorY()[fourthMDIndex] - mds.anchorY()[thirdMDIndex]));
1860  float sdOut_d = mds.anchorRt()[fourthMDIndex] - mds.anchorRt()[thirdMDIndex];
1861 
1862  runDeltaBetaIterationsT5(acc, betaIn, betaOut, betaAv, pt_beta, sdIn_dr, sdOut_dr, dr, lIn);
1863 
1864  const float betaInMMSF = (alpaka::math::abs(acc, betaInRHmin + betaInRHmax) > 0)
1865  ? (2.f * betaIn / alpaka::math::abs(acc, betaInRHmin + betaInRHmax))
1866  : 0.; //mean value of min,max is the old betaIn
1867  const float betaOutMMSF = (alpaka::math::abs(acc, betaOutRHmin + betaOutRHmax) > 0)
1868  ? (2.f * betaOut / alpaka::math::abs(acc, betaOutRHmin + betaOutRHmax))
1869  : 0.;
1870  betaInRHmin *= betaInMMSF;
1871  betaInRHmax *= betaInMMSF;
1872  betaOutRHmin *= betaOutMMSF;
1873  betaOutRHmax *= betaOutMMSF;
1874 
1875  float min_ptBeta_maxPtBeta = alpaka::math::min(
1876  acc, alpaka::math::abs(acc, pt_beta), kPt_betaMax); //need to confirm the range-out value of 7 GeV
1877  const float dBetaMuls2 = thetaMuls2 * 16.f / (min_ptBeta_maxPtBeta * min_ptBeta_maxPtBeta);
1878 
1879  const float alphaInAbsReg =
1880  alpaka::math::max(acc,
1881  alpaka::math::abs(acc, sdIn_alpha),
1882  alpaka::math::asin(acc, alpaka::math::min(acc, rt_InLo * k2Rinv1GeVf / 3.0f, kSinAlphaMax)));
1883  const float alphaOutAbsReg =
1884  alpaka::math::max(acc,
1885  alpaka::math::abs(acc, sdOut_alpha),
1886  alpaka::math::asin(acc, alpaka::math::min(acc, rt_OutLo * k2Rinv1GeVf / 3.0f, kSinAlphaMax)));
1887  const float dBetaInLum = lIn < 11 ? 0.0f : alpaka::math::abs(acc, alphaInAbsReg * kDeltaZLum / z_InLo);
1888  const float dBetaOutLum = lOut < 11 ? 0.0f : alpaka::math::abs(acc, alphaOutAbsReg * kDeltaZLum / z_OutLo);
1889  const float dBetaLum2 = (dBetaInLum + dBetaOutLum) * (dBetaInLum + dBetaOutLum);
1890 
1891  const float dBetaRIn2 = 0; // TODO-RH
1892 
1893  float dBetaROut2 = 0; //TODO-RH
1894  float betaOutCut = alpaka::math::asin(acc, alpaka::math::min(acc, dr * k2Rinv1GeVf / ptCut, kSinAlphaMax)) +
1895  (0.02f / sdOut_d) + alpaka::math::sqrt(acc, dBetaLum2 + dBetaMuls2);
1896 
1897  //Cut #6: The real beta cut
1898  if (alpaka::math::abs(acc, betaOut) >= betaOutCut)
1899  return false;
1900 
1901  float dBetaRes = 0.02f / alpaka::math::min(acc, sdOut_d, sdIn_d);
1902  float dBetaCut2 =
1903  (dBetaRes * dBetaRes * 2.0f + dBetaMuls2 + dBetaLum2 + dBetaRIn2 + dBetaROut2 +
1904  0.25f *
1905  (alpaka::math::abs(acc, betaInRHmin - betaInRHmax) + alpaka::math::abs(acc, betaOutRHmin - betaOutRHmax)) *
1906  (alpaka::math::abs(acc, betaInRHmin - betaInRHmax) + alpaka::math::abs(acc, betaOutRHmin - betaOutRHmax)));
1907  float dBeta = betaIn - betaOut;
1908  //Cut #7: Cut on dBeta
1909  return dBeta * dBeta <= dBetaCut2;
1910  }
1911 
1912  template <typename TAcc>
1913  ALPAKA_FN_ACC ALPAKA_FN_INLINE bool runQuintupletAlgoSelector(TAcc const& acc,
1914  ModulesConst modules,
1915  MiniDoubletsConst mds,
1916  SegmentsConst segments,
1917  uint16_t innerInnerLowerModuleIndex,
1918  uint16_t innerOuterLowerModuleIndex,
1919  uint16_t outerInnerLowerModuleIndex,
1920  uint16_t outerOuterLowerModuleIndex,
1921  unsigned int innerSegmentIndex,
1922  unsigned int outerSegmentIndex,
1923  unsigned int firstMDIndex,
1924  unsigned int secondMDIndex,
1925  unsigned int thirdMDIndex,
1926  unsigned int fourthMDIndex) {
1927  short innerInnerLowerModuleSubdet = modules.subdets()[innerInnerLowerModuleIndex];
1928  short innerOuterLowerModuleSubdet = modules.subdets()[innerOuterLowerModuleIndex];
1929  short outerInnerLowerModuleSubdet = modules.subdets()[outerInnerLowerModuleIndex];
1930  short outerOuterLowerModuleSubdet = modules.subdets()[outerOuterLowerModuleIndex];
1931 
1932  if (innerInnerLowerModuleSubdet == Barrel and innerOuterLowerModuleSubdet == Barrel and
1933  outerInnerLowerModuleSubdet == Barrel and outerOuterLowerModuleSubdet == Barrel) {
1934  return runQuintupletDefaultAlgoBBBB(acc,
1935  modules,
1936  mds,
1937  segments,
1938  innerInnerLowerModuleIndex,
1939  innerOuterLowerModuleIndex,
1940  outerInnerLowerModuleIndex,
1941  outerOuterLowerModuleIndex,
1942  innerSegmentIndex,
1943  outerSegmentIndex,
1944  firstMDIndex,
1945  secondMDIndex,
1946  thirdMDIndex,
1947  fourthMDIndex);
1948  } else if (innerInnerLowerModuleSubdet == Barrel and innerOuterLowerModuleSubdet == Barrel and
1949  outerInnerLowerModuleSubdet == Endcap and outerOuterLowerModuleSubdet == Endcap) {
1950  return runQuintupletDefaultAlgoBBEE(acc,
1951  modules,
1952  mds,
1953  segments,
1954  innerInnerLowerModuleIndex,
1955  innerOuterLowerModuleIndex,
1956  outerInnerLowerModuleIndex,
1957  outerOuterLowerModuleIndex,
1958  innerSegmentIndex,
1959  outerSegmentIndex,
1960  firstMDIndex,
1961  secondMDIndex,
1962  thirdMDIndex,
1963  fourthMDIndex);
1964  } else if (innerInnerLowerModuleSubdet == Barrel and innerOuterLowerModuleSubdet == Barrel and
1965  outerInnerLowerModuleSubdet == Barrel and outerOuterLowerModuleSubdet == Endcap) {
1966  return runQuintupletDefaultAlgoBBBB(acc,
1967  modules,
1968  mds,
1969  segments,
1970  innerInnerLowerModuleIndex,
1971  innerOuterLowerModuleIndex,
1972  outerInnerLowerModuleIndex,
1973  outerOuterLowerModuleIndex,
1974  innerSegmentIndex,
1975  outerSegmentIndex,
1976  firstMDIndex,
1977  secondMDIndex,
1978  thirdMDIndex,
1979  fourthMDIndex);
1980  } else if (innerInnerLowerModuleSubdet == Barrel and innerOuterLowerModuleSubdet == Endcap and
1981  outerInnerLowerModuleSubdet == Endcap and outerOuterLowerModuleSubdet == Endcap) {
1982  return runQuintupletDefaultAlgoBBEE(acc,
1983  modules,
1984  mds,
1985  segments,
1986  innerInnerLowerModuleIndex,
1987  innerOuterLowerModuleIndex,
1988  outerInnerLowerModuleIndex,
1989  outerOuterLowerModuleIndex,
1990  innerSegmentIndex,
1991  outerSegmentIndex,
1992  firstMDIndex,
1993  secondMDIndex,
1994  thirdMDIndex,
1995  fourthMDIndex);
1996  } else if (innerInnerLowerModuleSubdet == Endcap and innerOuterLowerModuleSubdet == Endcap and
1997  outerInnerLowerModuleSubdet == Endcap and outerOuterLowerModuleSubdet == Endcap) {
1998  return runQuintupletDefaultAlgoEEEE(acc,
1999  modules,
2000  mds,
2001  segments,
2002  innerInnerLowerModuleIndex,
2003  innerOuterLowerModuleIndex,
2004  outerInnerLowerModuleIndex,
2005  outerOuterLowerModuleIndex,
2006  innerSegmentIndex,
2007  outerSegmentIndex,
2008  firstMDIndex,
2009  secondMDIndex,
2010  thirdMDIndex,
2011  fourthMDIndex);
2012  }
2013 
2014  return false;
2015  }
2016 
2017  template <typename TAcc>
2018  ALPAKA_FN_ACC ALPAKA_FN_INLINE bool runQuintupletDefaultAlgo(TAcc const& acc,
2019  ModulesConst modules,
2020  MiniDoubletsConst mds,
2021  SegmentsConst segments,
2022  TripletsConst triplets,
2023  uint16_t lowerModuleIndex1,
2024  uint16_t lowerModuleIndex2,
2025  uint16_t lowerModuleIndex3,
2026  uint16_t lowerModuleIndex4,
2027  uint16_t lowerModuleIndex5,
2028  unsigned int innerTripletIndex,
2029  unsigned int outerTripletIndex,
2030  float& innerRadius,
2031  float& outerRadius,
2032  float& bridgeRadius,
2033  float& regressionG,
2034  float& regressionF,
2035  float& regressionRadius,
2036  float& rzChiSquared,
2037  float& chiSquared,
2038  float& nonAnchorChiSquared,
2039  bool& tightCutFlag) {
2040  unsigned int firstSegmentIndex = triplets.segmentIndices()[innerTripletIndex][0];
2041  unsigned int secondSegmentIndex = triplets.segmentIndices()[innerTripletIndex][1];
2042  unsigned int thirdSegmentIndex = triplets.segmentIndices()[outerTripletIndex][0];
2043  unsigned int fourthSegmentIndex = triplets.segmentIndices()[outerTripletIndex][1];
2044 
2045  unsigned int innerOuterOuterMiniDoubletIndex =
2046  segments.mdIndices()[secondSegmentIndex][1]; //inner triplet outer segment outer MD index
2047  unsigned int outerInnerInnerMiniDoubletIndex =
2048  segments.mdIndices()[thirdSegmentIndex][0]; //outer triplet inner segment inner MD index
2049 
2050  //this cut reduces the number of candidates by a factor of 3, i.e., 2 out of 3 warps can end right here!
2051  if (innerOuterOuterMiniDoubletIndex != outerInnerInnerMiniDoubletIndex)
2052  return false;
2053 
2054  unsigned int firstMDIndex = segments.mdIndices()[firstSegmentIndex][0];
2055  unsigned int secondMDIndex = segments.mdIndices()[secondSegmentIndex][0];
2056  unsigned int thirdMDIndex = segments.mdIndices()[secondSegmentIndex][1];
2057  unsigned int fourthMDIndex = segments.mdIndices()[thirdSegmentIndex][1];
2058  unsigned int fifthMDIndex = segments.mdIndices()[fourthSegmentIndex][1];
2059 
2060  if (not runQuintupletAlgoSelector(acc,
2061  modules,
2062  mds,
2063  segments,
2064  lowerModuleIndex1,
2065  lowerModuleIndex2,
2066  lowerModuleIndex3,
2067  lowerModuleIndex4,
2068  firstSegmentIndex,
2069  thirdSegmentIndex,
2070  firstMDIndex,
2071  secondMDIndex,
2072  thirdMDIndex,
2073  fourthMDIndex))
2074  return false;
2075 
2076  if (not runQuintupletAlgoSelector(acc,
2077  modules,
2078  mds,
2079  segments,
2080  lowerModuleIndex1,
2081  lowerModuleIndex2,
2082  lowerModuleIndex4,
2083  lowerModuleIndex5,
2084  firstSegmentIndex,
2085  fourthSegmentIndex,
2086  firstMDIndex,
2087  secondMDIndex,
2088  fourthMDIndex,
2089  fifthMDIndex))
2090  return false;
2091 
2092  float x1 = mds.anchorX()[firstMDIndex];
2093  float x2 = mds.anchorX()[secondMDIndex];
2094  float x3 = mds.anchorX()[thirdMDIndex];
2095  float x4 = mds.anchorX()[fourthMDIndex];
2096  float x5 = mds.anchorX()[fifthMDIndex];
2097 
2098  float y1 = mds.anchorY()[firstMDIndex];
2099  float y2 = mds.anchorY()[secondMDIndex];
2100  float y3 = mds.anchorY()[thirdMDIndex];
2101  float y4 = mds.anchorY()[fourthMDIndex];
2102  float y5 = mds.anchorY()[fifthMDIndex];
2103 
2104  //construct the arrays
2105  float x1Vec[] = {x1, x1, x1};
2106  float y1Vec[] = {y1, y1, y1};
2107  float x2Vec[] = {x2, x2, x2};
2108  float y2Vec[] = {y2, y2, y2};
2109  float x3Vec[] = {x3, x3, x3};
2110  float y3Vec[] = {y3, y3, y3};
2111 
2112  if (modules.subdets()[lowerModuleIndex1] == Endcap and modules.moduleType()[lowerModuleIndex1] == TwoS) {
2113  x1Vec[1] = mds.anchorLowEdgeX()[firstMDIndex];
2114  x1Vec[2] = mds.anchorHighEdgeX()[firstMDIndex];
2115 
2116  y1Vec[1] = mds.anchorLowEdgeY()[firstMDIndex];
2117  y1Vec[2] = mds.anchorHighEdgeY()[firstMDIndex];
2118  }
2119  if (modules.subdets()[lowerModuleIndex2] == Endcap and modules.moduleType()[lowerModuleIndex2] == TwoS) {
2120  x2Vec[1] = mds.anchorLowEdgeX()[secondMDIndex];
2121  x2Vec[2] = mds.anchorHighEdgeX()[secondMDIndex];
2122 
2123  y2Vec[1] = mds.anchorLowEdgeY()[secondMDIndex];
2124  y2Vec[2] = mds.anchorHighEdgeY()[secondMDIndex];
2125  }
2126  if (modules.subdets()[lowerModuleIndex3] == Endcap and modules.moduleType()[lowerModuleIndex3] == TwoS) {
2127  x3Vec[1] = mds.anchorLowEdgeX()[thirdMDIndex];
2128  x3Vec[2] = mds.anchorHighEdgeX()[thirdMDIndex];
2129 
2130  y3Vec[1] = mds.anchorLowEdgeY()[thirdMDIndex];
2131  y3Vec[2] = mds.anchorHighEdgeY()[thirdMDIndex];
2132  }
2133 
2134  float innerRadiusMin2S, innerRadiusMax2S;
2135  computeErrorInRadius(acc, x1Vec, y1Vec, x2Vec, y2Vec, x3Vec, y3Vec, innerRadiusMin2S, innerRadiusMax2S);
2136 
2137  for (int i = 0; i < 3; i++) {
2138  x1Vec[i] = x4;
2139  y1Vec[i] = y4;
2140  }
2141  if (modules.subdets()[lowerModuleIndex4] == Endcap and modules.moduleType()[lowerModuleIndex4] == TwoS) {
2142  x1Vec[1] = mds.anchorLowEdgeX()[fourthMDIndex];
2143  x1Vec[2] = mds.anchorHighEdgeX()[fourthMDIndex];
2144 
2145  y1Vec[1] = mds.anchorLowEdgeY()[fourthMDIndex];
2146  y1Vec[2] = mds.anchorHighEdgeY()[fourthMDIndex];
2147  }
2148 
2149  float bridgeRadiusMin2S, bridgeRadiusMax2S;
2150  computeErrorInRadius(acc, x2Vec, y2Vec, x3Vec, y3Vec, x1Vec, y1Vec, bridgeRadiusMin2S, bridgeRadiusMax2S);
2151 
2152  for (int i = 0; i < 3; i++) {
2153  x2Vec[i] = x5;
2154  y2Vec[i] = y5;
2155  }
2156  if (modules.subdets()[lowerModuleIndex5] == Endcap and modules.moduleType()[lowerModuleIndex5] == TwoS) {
2157  x2Vec[1] = mds.anchorLowEdgeX()[fifthMDIndex];
2158  x2Vec[2] = mds.anchorHighEdgeX()[fifthMDIndex];
2159 
2160  y2Vec[1] = mds.anchorLowEdgeY()[fifthMDIndex];
2161  y2Vec[2] = mds.anchorHighEdgeY()[fifthMDIndex];
2162  }
2163 
2164  float outerRadiusMin2S, outerRadiusMax2S;
2165  computeErrorInRadius(acc, x3Vec, y3Vec, x1Vec, y1Vec, x2Vec, y2Vec, outerRadiusMin2S, outerRadiusMax2S);
2166 
2167  float g, f;
2168  outerRadius = triplets.radius()[outerTripletIndex];
2169  bridgeRadius = computeRadiusFromThreeAnchorHits(acc, x2, y2, x3, y3, x4, y4, g, f);
2170  innerRadius = triplets.radius()[innerTripletIndex];
2171  g = triplets.centerX()[innerTripletIndex];
2172  f = triplets.centerY()[innerTripletIndex];
2173 
2174  float inner_pt = 2 * k2Rinv1GeVf * innerRadius;
2175 
2176  if (not passT5RZConstraint(acc,
2177  modules,
2178  mds,
2179  firstMDIndex,
2180  secondMDIndex,
2181  thirdMDIndex,
2182  fourthMDIndex,
2183  fifthMDIndex,
2184  lowerModuleIndex1,
2185  lowerModuleIndex2,
2186  lowerModuleIndex3,
2187  lowerModuleIndex4,
2188  lowerModuleIndex5,
2189  rzChiSquared,
2190  inner_pt,
2191  innerRadius,
2192  g,
2193  f,
2194  tightCutFlag))
2195  return false;
2196 
2197  if (innerRadius < 0.95f * ptCut / (2.f * k2Rinv1GeVf))
2198  return false;
2199 
2200  //split by category
2201  bool matchedRadii;
2202  if (modules.subdets()[lowerModuleIndex1] == Barrel and modules.subdets()[lowerModuleIndex2] == Barrel and
2203  modules.subdets()[lowerModuleIndex3] == Barrel and modules.subdets()[lowerModuleIndex4] == Barrel and
2204  modules.subdets()[lowerModuleIndex5] == Barrel) {
2205  matchedRadii = matchRadiiBBBBB(acc, innerRadius, bridgeRadius, outerRadius);
2206  } else if (modules.subdets()[lowerModuleIndex1] == Barrel and modules.subdets()[lowerModuleIndex2] == Barrel and
2207  modules.subdets()[lowerModuleIndex3] == Barrel and modules.subdets()[lowerModuleIndex4] == Barrel and
2208  modules.subdets()[lowerModuleIndex5] == Endcap) {
2209  matchedRadii = matchRadiiBBBBE(acc, innerRadius, bridgeRadius, outerRadius);
2210  } else if (modules.subdets()[lowerModuleIndex1] == Barrel and modules.subdets()[lowerModuleIndex2] == Barrel and
2211  modules.subdets()[lowerModuleIndex3] == Barrel and modules.subdets()[lowerModuleIndex4] == Endcap and
2212  modules.subdets()[lowerModuleIndex5] == Endcap) {
2213  if (modules.layers()[lowerModuleIndex1] == 1) {
2214  matchedRadii =
2215  matchRadiiBBBEE12378(acc, innerRadius, bridgeRadius, outerRadius, bridgeRadiusMin2S, bridgeRadiusMax2S);
2216  } else if (modules.layers()[lowerModuleIndex1] == 2) {
2217  matchedRadii =
2218  matchRadiiBBBEE23478(acc, innerRadius, bridgeRadius, outerRadius, bridgeRadiusMin2S, bridgeRadiusMax2S);
2219  } else {
2220  matchedRadii =
2221  matchRadiiBBBEE34578(acc, innerRadius, bridgeRadius, outerRadius, bridgeRadiusMin2S, bridgeRadiusMax2S);
2222  }
2223  }
2224 
2225  else if (modules.subdets()[lowerModuleIndex1] == Barrel and modules.subdets()[lowerModuleIndex2] == Barrel and
2226  modules.subdets()[lowerModuleIndex3] == Endcap and modules.subdets()[lowerModuleIndex4] == Endcap and
2227  modules.subdets()[lowerModuleIndex5] == Endcap) {
2228  matchedRadii = matchRadiiBBEEE(acc, innerRadius, bridgeRadius, outerRadius, bridgeRadiusMin2S, bridgeRadiusMax2S);
2229  } else if (modules.subdets()[lowerModuleIndex1] == Barrel and modules.subdets()[lowerModuleIndex2] == Endcap and
2230  modules.subdets()[lowerModuleIndex3] == Endcap and modules.subdets()[lowerModuleIndex4] == Endcap and
2231  modules.subdets()[lowerModuleIndex5] == Endcap) {
2232  matchedRadii = matchRadiiBEEEE(acc,
2233  innerRadius,
2234  bridgeRadius,
2235  outerRadius,
2236  innerRadiusMin2S,
2237  innerRadiusMax2S,
2238  bridgeRadiusMin2S,
2239  bridgeRadiusMax2S);
2240  } else {
2241  matchedRadii = matchRadiiEEEEE(acc,
2242  innerRadius,
2243  bridgeRadius,
2244  outerRadius,
2245  innerRadiusMin2S,
2246  innerRadiusMax2S,
2247  bridgeRadiusMin2S,
2248  bridgeRadiusMax2S);
2249  }
2250 
2251  //compute regression radius right here - this computation is expensive!!!
2252  if (not matchedRadii)
2253  return false;
2254 
2255  float xVec[] = {x1, x2, x3, x4, x5};
2256  float yVec[] = {y1, y2, y3, y4, y5};
2257  const uint16_t lowerModuleIndices[] = {
2258  lowerModuleIndex1, lowerModuleIndex2, lowerModuleIndex3, lowerModuleIndex4, lowerModuleIndex5};
2259 
2260  // 5 categories for sigmas
2261  float sigmas2[5], delta1[5], delta2[5], slopes[5];
2262  bool isFlat[5];
2263 
2264  computeSigmasForRegression(acc, modules, lowerModuleIndices, delta1, delta2, slopes, isFlat);
2265  regressionRadius = computeRadiusUsingRegression(acc,
2266  Params_T5::kLayers,
2267  xVec,
2268  yVec,
2269  delta1,
2270  delta2,
2271  slopes,
2272  isFlat,
2273  regressionG,
2274  regressionF,
2275  sigmas2,
2276  chiSquared);
2277 
2278  unsigned int mdIndices[] = {firstMDIndex, secondMDIndex, thirdMDIndex, fourthMDIndex, fifthMDIndex};
2279  float inference = t5dnn::runInference(acc,
2280  modules,
2281  mds,
2282  segments,
2283  triplets,
2284  xVec,
2285  yVec,
2286  mdIndices,
2287  lowerModuleIndices,
2288  innerTripletIndex,
2289  outerTripletIndex,
2290  innerRadius,
2291  outerRadius,
2292  bridgeRadius);
2293  tightCutFlag = tightCutFlag and (inference > t5dnn::kLSTWp2); // T5-in-TC cut
2294  if (inference <= t5dnn::kLSTWp2) // T5-building cut
2295  return false;
2296 
2297  //compute the other chisquared
2298  //non anchor is always shifted for tilted and endcap!
2299  float nonAnchorDelta1[Params_T5::kLayers], nonAnchorDelta2[Params_T5::kLayers], nonAnchorSlopes[Params_T5::kLayers];
2300  float nonAnchorxs[] = {mds.outerX()[firstMDIndex],
2301  mds.outerX()[secondMDIndex],
2302  mds.outerX()[thirdMDIndex],
2303  mds.outerX()[fourthMDIndex],
2304  mds.outerX()[fifthMDIndex]};
2305  float nonAnchorys[] = {mds.outerY()[firstMDIndex],
2306  mds.outerY()[secondMDIndex],
2307  mds.outerY()[thirdMDIndex],
2308  mds.outerY()[fourthMDIndex],
2309  mds.outerY()[fifthMDIndex]};
2310 
2311  computeSigmasForRegression(acc,
2312  modules,
2313  lowerModuleIndices,
2314  nonAnchorDelta1,
2315  nonAnchorDelta2,
2316  nonAnchorSlopes,
2317  isFlat,
2318  Params_T5::kLayers,
2319  false);
2320  nonAnchorChiSquared = computeChiSquared(acc,
2321  Params_T5::kLayers,
2322  nonAnchorxs,
2323  nonAnchorys,
2324  nonAnchorDelta1,
2325  nonAnchorDelta2,
2326  nonAnchorSlopes,
2327  isFlat,
2328  regressionG,
2329  regressionF,
2330  regressionRadius);
2331  return true;
2332  }
2333 
2334  struct CreateQuintuplets {
2335  template <typename TAcc>
2336  ALPAKA_FN_ACC void operator()(TAcc const& acc,
2337  ModulesConst modules,
2338  MiniDoubletsConst mds,
2339  SegmentsConst segments,
2340  Triplets triplets,
2341  TripletsOccupancyConst tripletsOccupancy,
2342  Quintuplets quintuplets,
2343  QuintupletsOccupancy quintupletsOccupancy,
2344  ObjectRangesConst ranges,
2345  uint16_t nEligibleT5Modules) const {
2346  auto const globalThreadIdx = alpaka::getIdx<alpaka::Grid, alpaka::Threads>(acc);
2347  auto const gridThreadExtent = alpaka::getWorkDiv<alpaka::Grid, alpaka::Threads>(acc);
2348 
2349  for (int iter = globalThreadIdx[0]; iter < nEligibleT5Modules; iter += gridThreadExtent[0]) {
2350  uint16_t lowerModule1 = ranges.indicesOfEligibleT5Modules()[iter];
2351  short layer2_adjustment;
2352  int layer = modules.layers()[lowerModule1];
2353  if (layer == 1) {
2354  layer2_adjustment = 1;
2355  } // get upper segment to be in second layer
2356  else if (layer == 2) {
2357  layer2_adjustment = 0;
2358  } // get lower segment to be in second layer
2359  else {
2360  continue;
2361  }
2362  unsigned int nInnerTriplets = tripletsOccupancy.nTriplets()[lowerModule1];
2363  for (unsigned int innerTripletArrayIndex = globalThreadIdx[1]; innerTripletArrayIndex < nInnerTriplets;
2364  innerTripletArrayIndex += gridThreadExtent[1]) {
2365  unsigned int innerTripletIndex = ranges.tripletModuleIndices()[lowerModule1] + innerTripletArrayIndex;
2366  uint16_t lowerModule2 = triplets.lowerModuleIndices()[innerTripletIndex][1];
2367  uint16_t lowerModule3 = triplets.lowerModuleIndices()[innerTripletIndex][2];
2368  unsigned int nOuterTriplets = tripletsOccupancy.nTriplets()[lowerModule3];
2369  for (unsigned int outerTripletArrayIndex = globalThreadIdx[2]; outerTripletArrayIndex < nOuterTriplets;
2370  outerTripletArrayIndex += gridThreadExtent[2]) {
2371  unsigned int outerTripletIndex = ranges.tripletModuleIndices()[lowerModule3] + outerTripletArrayIndex;
2372  uint16_t lowerModule4 = triplets.lowerModuleIndices()[outerTripletIndex][1];
2373  uint16_t lowerModule5 = triplets.lowerModuleIndices()[outerTripletIndex][2];
2374 
2375  float innerRadius, outerRadius, bridgeRadius, regressionG, regressionF, regressionRadius, rzChiSquared,
2376  chiSquared, nonAnchorChiSquared; //required for making distributions
2377 
2378  bool tightCutFlag = false;
2379  bool success = runQuintupletDefaultAlgo(acc,
2380  modules,
2381  mds,
2382  segments,
2383  triplets,
2384  lowerModule1,
2385  lowerModule2,
2386  lowerModule3,
2387  lowerModule4,
2388  lowerModule5,
2389  innerTripletIndex,
2390  outerTripletIndex,
2391  innerRadius,
2392  outerRadius,
2393  bridgeRadius,
2394  regressionG,
2395  regressionF,
2396  regressionRadius,
2397  rzChiSquared,
2398  chiSquared,
2399  nonAnchorChiSquared,
2400  tightCutFlag);
2401 
2402  if (success) {
2403  int totOccupancyQuintuplets = alpaka::atomicAdd(
2404  acc, &quintupletsOccupancy.totOccupancyQuintuplets()[lowerModule1], 1u, alpaka::hierarchy::Threads{});
2405  if (totOccupancyQuintuplets >= ranges.quintupletModuleOccupancy()[lowerModule1]) {
2406 #ifdef WARNINGS
2407  printf("Quintuplet excess alert! Module index = %d\n", lowerModule1);
2408 #endif
2409  } else {
2410  int quintupletModuleIndex = alpaka::atomicAdd(
2411  acc, &quintupletsOccupancy.nQuintuplets()[lowerModule1], 1u, alpaka::hierarchy::Threads{});
2412  //this if statement should never get executed!
2413  if (ranges.quintupletModuleIndices()[lowerModule1] == -1) {
2414 #ifdef WARNINGS
2415  printf("Quintuplets : no memory for module at module index = %d\n", lowerModule1);
2416 #endif
2417  } else {
2418  unsigned int quintupletIndex = ranges.quintupletModuleIndices()[lowerModule1] + quintupletModuleIndex;
2419  float phi = mds.anchorPhi()[segments.mdIndices()[triplets.segmentIndices()[innerTripletIndex][0]]
2420  [layer2_adjustment]];
2421  float eta = mds.anchorEta()[segments.mdIndices()[triplets.segmentIndices()[innerTripletIndex][0]]
2422  [layer2_adjustment]];
2423  float pt = (innerRadius + outerRadius) * k2Rinv1GeVf;
2424  float scores = chiSquared + nonAnchorChiSquared;
2425  addQuintupletToMemory(triplets,
2426  quintuplets,
2427  innerTripletIndex,
2428  outerTripletIndex,
2429  lowerModule1,
2430  lowerModule2,
2431  lowerModule3,
2432  lowerModule4,
2433  lowerModule5,
2434  innerRadius,
2435  bridgeRadius,
2436  outerRadius,
2437  regressionG,
2438  regressionF,
2439  regressionRadius,
2440  rzChiSquared,
2441  chiSquared,
2442  nonAnchorChiSquared,
2443  pt,
2444  eta,
2445  phi,
2446  scores,
2447  layer,
2448  quintupletIndex,
2449  tightCutFlag);
2450 
2451  triplets.partOfT5()[quintuplets.tripletIndices()[quintupletIndex][0]] = true;
2452  triplets.partOfT5()[quintuplets.tripletIndices()[quintupletIndex][1]] = true;
2453  }
2454  }
2455  }
2456  }
2457  }
2458  }
2459  }
2460  };
2461 
2462  struct CreateEligibleModulesListForQuintuplets {
2463  template <typename TAcc>
2464  ALPAKA_FN_ACC void operator()(TAcc const& acc,
2465  ModulesConst modules,
2466  TripletsOccupancyConst tripletsOccupancy,
2467  ObjectRanges ranges) const {
2468  // implementation is 1D with a single block
2469  static_assert(std::is_same_v<TAcc, ALPAKA_ACCELERATOR_NAMESPACE::Acc1D>, "Should be Acc1D");
2470  ALPAKA_ASSERT_ACC((alpaka::getWorkDiv<alpaka::Grid, alpaka::Blocks>(acc)[0] == 1));
2471 
2472  auto const globalThreadIdx = alpaka::getIdx<alpaka::Grid, alpaka::Threads>(acc);
2473  auto const gridThreadExtent = alpaka::getWorkDiv<alpaka::Grid, alpaka::Threads>(acc);
2474 
2475  // Initialize variables in shared memory and set to 0
2476  int& nEligibleT5Modulesx = alpaka::declareSharedVar<int, __COUNTER__>(acc);
2477  int& nTotalQuintupletsx = alpaka::declareSharedVar<int, __COUNTER__>(acc);
2478  if (cms::alpakatools::once_per_block(acc)) {
2479  nTotalQuintupletsx = 0;
2480  nEligibleT5Modulesx = 0;
2481  }
2482  alpaka::syncBlockThreads(acc);
2483 
2484  for (int i = globalThreadIdx[0]; i < modules.nLowerModules(); i += gridThreadExtent[0]) {
2485  // Condition for a quintuple to exist for a module
2486  // TCs don't exist for layers 5 and 6 barrel, and layers 2,3,4,5 endcap
2487  short module_rings = modules.rings()[i];
2488  short module_layers = modules.layers()[i];
2489  short module_subdets = modules.subdets()[i];
2490  float module_eta = alpaka::math::abs(acc, modules.eta()[i]);
2491 
2492  if (tripletsOccupancy.nTriplets()[i] == 0)
2493  continue;
2494  if (module_subdets == Barrel and module_layers >= 3)
2495  continue;
2496  if (module_subdets == Endcap and module_layers > 1)
2497  continue;
2498 
2499  int nEligibleT5Modules = alpaka::atomicAdd(acc, &nEligibleT5Modulesx, 1, alpaka::hierarchy::Threads{});
2500 
2501  int category_number;
2502  if (module_layers <= 3 && module_subdets == 5)
2503  category_number = 0;
2504  else if (module_layers >= 4 && module_subdets == 5)
2505  category_number = 1;
2506  else if (module_layers <= 2 && module_subdets == 4 && module_rings >= 11)
2507  category_number = 2;
2508  else if (module_layers >= 3 && module_subdets == 4 && module_rings >= 8)
2509  category_number = 2;
2510  else if (module_layers <= 2 && module_subdets == 4 && module_rings <= 10)
2511  category_number = 3;
2512  else if (module_layers >= 3 && module_subdets == 4 && module_rings <= 7)
2513  category_number = 3;
2514  else
2515  category_number = -1;
2516 
2517  int eta_number;
2518  if (module_eta < 0.75f)
2519  eta_number = 0;
2520  else if (module_eta < 1.5f)
2521  eta_number = 1;
2522  else if (module_eta < 2.25f)
2523  eta_number = 2;
2524  else if (module_eta < 3.0f)
2525  eta_number = 3;
2526  else
2527  eta_number = -1;
2528 
2529  int occupancy;
2530  if (category_number == 0 && eta_number == 0)
2531  occupancy = 336;
2532  else if (category_number == 0 && eta_number == 1)
2533  occupancy = 414;
2534  else if (category_number == 0 && eta_number == 2)
2535  occupancy = 231;
2536  else if (category_number == 0 && eta_number == 3)
2537  occupancy = 146;
2538  else if (category_number == 3 && eta_number == 1)
2539  occupancy = 0;
2540  else if (category_number == 3 && eta_number == 2)
2541  occupancy = 191;
2542  else if (category_number == 3 && eta_number == 3)
2543  occupancy = 106;
2544  else {
2545  occupancy = 0;
2546 #ifdef WARNINGS
2547  printf("Unhandled case in createEligibleModulesListForQuintupletsGPU! Module index = %i\n", i);
2548 #endif
2549  }
2550 
2551  int nTotQ = alpaka::atomicAdd(acc, &nTotalQuintupletsx, occupancy, alpaka::hierarchy::Threads{});
2552  ranges.quintupletModuleIndices()[i] = nTotQ;
2553  ranges.indicesOfEligibleT5Modules()[nEligibleT5Modules] = i;
2554  ranges.quintupletModuleOccupancy()[i] = occupancy;
2555  }
2556 
2557  // Wait for all threads to finish before reporting final values
2558  alpaka::syncBlockThreads(acc);
2559  if (cms::alpakatools::once_per_block(acc)) {
2560  ranges.nEligibleT5Modules() = static_cast<uint16_t>(nEligibleT5Modulesx);
2561  ranges.nTotalQuints() = static_cast<unsigned int>(nTotalQuintupletsx);
2562  }
2563  }
2564  };
2565 
2566  struct AddQuintupletRangesToEventExplicit {
2567  template <typename TAcc>
2568  ALPAKA_FN_ACC void operator()(TAcc const& acc,
2569  ModulesConst modules,
2570  QuintupletsOccupancyConst quintupletsOccupancy,
2571  ObjectRanges ranges) const {
2572  // implementation is 1D with a single block
2573  static_assert(std::is_same_v<TAcc, ALPAKA_ACCELERATOR_NAMESPACE::Acc1D>, "Should be Acc1D");
2574  ALPAKA_ASSERT_ACC((alpaka::getWorkDiv<alpaka::Grid, alpaka::Blocks>(acc)[0] == 1));
2575 
2576  auto const globalThreadIdx = alpaka::getIdx<alpaka::Grid, alpaka::Threads>(acc);
2577  auto const gridThreadExtent = alpaka::getWorkDiv<alpaka::Grid, alpaka::Threads>(acc);
2578 
2579  for (uint16_t i = globalThreadIdx[0]; i < modules.nLowerModules(); i += gridThreadExtent[0]) {
2580  if (quintupletsOccupancy.nQuintuplets()[i] == 0 or ranges.quintupletModuleIndices()[i] == -1) {
2581  ranges.quintupletRanges()[i][0] = -1;
2582  ranges.quintupletRanges()[i][1] = -1;
2583  } else {
2584  ranges.quintupletRanges()[i][0] = ranges.quintupletModuleIndices()[i];
2585  ranges.quintupletRanges()[i][1] =
2586  ranges.quintupletModuleIndices()[i] + quintupletsOccupancy.nQuintuplets()[i] - 1;
2587  }
2588  }
2589  }
2590  };
2591 } // namespace ALPAKA_ACCELERATOR_NAMESPACE::lst
2592 #endif
ppsModifySingularModes_cfi.z2
z2
Definition: ppsModifySingularModes_cfi.py:9
ALPAKA_ACCELERATOR_NAMESPACE::lst::t5dnn::runInference
ALPAKA_FN_ACC ALPAKA_FN_INLINE float runInference(TAcc const &acc, ModulesConst modules, MiniDoubletsConst mds, SegmentsConst segments, TripletsConst triplets, const float *xVec, const float *yVec, const unsigned int *mdIndices, const uint16_t *lowerModuleIndices, unsigned int innerTripletIndex, unsigned int outerTripletIndex, float innerRadius, float outerRadius, float bridgeRadius)
Definition: NeuralNetwork.h:18
ALPAKA_ACCELERATOR_NAMESPACE::lst::runQuintupletAlgoSelector
ALPAKA_FN_ACC ALPAKA_FN_INLINE bool runQuintupletAlgoSelector(TAcc const &acc, ModulesConst modules, MiniDoubletsConst mds, SegmentsConst segments, uint16_t innerInnerLowerModuleIndex, uint16_t innerOuterLowerModuleIndex, uint16_t outerInnerLowerModuleIndex, uint16_t outerOuterLowerModuleIndex, unsigned int innerSegmentIndex, unsigned int outerSegmentIndex, unsigned int firstMDIndex, unsigned int secondMDIndex, unsigned int thirdMDIndex, unsigned int fourthMDIndex)
Definition: Quintuplet.h:1913
CommonMethods.isnan
def isnan(num)
Definition: CommonMethods.py:97
cms::alpakatools::once_per_block
ALPAKA_FN_ACC constexpr bool once_per_block(TAcc const &acc)
Definition: workdivision.h:1394
mps_fire.i
i
Definition: mps_fire.py:429
ALPAKA_ACCELERATOR_NAMESPACE::lst
Definition: Common.h:8
ALPAKA_ACCELERATOR_NAMESPACE::lst::passChiSquaredConstraint
ALPAKA_FN_ACC ALPAKA_FN_INLINE bool passChiSquaredConstraint(ModulesConst modules, uint16_t lowerModuleIndex1, uint16_t lowerModuleIndex2, uint16_t lowerModuleIndex3, uint16_t lowerModuleIndex4, uint16_t lowerModuleIndex5, float chiSquared)
Definition: Quintuplet.h:95
B
Definition: APVGainStruct.h:7
ALPAKA_ACCELERATOR_NAMESPACE::lst::computeSigmasForRegression
ALPAKA_FN_ACC ALPAKA_FN_INLINE void computeSigmasForRegression(TAcc const &acc, ModulesConst modules, const uint16_t *lowerModuleIndices, float *delta1, float *delta2, float *slopes, bool *isFlat, unsigned int nPoints=5, bool anchorHits=true)
Definition: Quintuplet.h:871
ALPAKA_ACCELERATOR_NAMESPACE::lst::checkIntervalOverlap
ALPAKA_FN_ACC ALPAKA_FN_INLINE bool checkIntervalOverlap(float firstMin, float firstMax, float secondMin, float secondMax)
Definition: Quintuplet.h:21
ALPAKA_ACCELERATOR_NAMESPACE::lst::CreateQuintuplets
Definition: Quintuplet.h:2334
ALPAKA_ACCELERATOR_NAMESPACE::lst::k2Rinv1GeVf
ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float k2Rinv1GeVf
Definition: Common.h:49
lst::ObjectRanges
ObjectRangesSoA::View ObjectRanges
Definition: ObjectRangesSoA.h:31
WZElectronSkims53X_cff.max
max
Definition: WZElectronSkims53X_cff.py:197
slope
static const double slope[3]
Definition: CastorTimeSlew.cc:6
correctedFixedMatrixSuperClustersWithPreshower_cfi.corrF
corrF
Definition: correctedFixedMatrixSuperClustersWithPreshower_cfi.py:19
tas::module_rings
const std::vector< int > & module_rings()
Definition: LSTEff.cc:3384
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
ALPAKA_ACCELERATOR_NAMESPACE::lst::runQuintupletDefaultAlgo
ALPAKA_FN_ACC ALPAKA_FN_INLINE bool runQuintupletDefaultAlgo(TAcc const &acc, ModulesConst modules, MiniDoubletsConst mds, SegmentsConst segments, TripletsConst triplets, uint16_t lowerModuleIndex1, uint16_t lowerModuleIndex2, uint16_t lowerModuleIndex3, uint16_t lowerModuleIndex4, uint16_t lowerModuleIndex5, unsigned int innerTripletIndex, unsigned int outerTripletIndex, float &innerRadius, float &outerRadius, float &bridgeRadius, float &regressionG, float &regressionF, float &regressionRadius, float &rzChiSquared, float &chiSquared, float &nonAnchorChiSquared, bool &tightCutFlag)
Definition: Quintuplet.h:2018
testProducerWithPsetDescEmpty_cfi.y2
y2
Definition: testProducerWithPsetDescEmpty_cfi.py:30
ALPAKA_ACCELERATOR_NAMESPACE::ecal::reconstruction::internal::barrel::side
ALPAKA_FN_ACC int side(int ieta, int iphi)
Definition: KernelHelpers.dev.cc:92
ALPAKA_ACCELERATOR_NAMESPACE::lst::kWidthPS
ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kWidthPS
Definition: Common.h:58
ALPAKA_ACCELERATOR_NAMESPACE::lst::computeRadiusFromThreeAnchorHits
ALPAKA_FN_ACC ALPAKA_FN_INLINE float computeRadiusFromThreeAnchorHits(TAcc const &acc, float x1, float y1, float x2, float y2, float x3, float y3, float &g, float &f)
Definition: Triplet.h:567
ALPAKA_ACCELERATOR_NAMESPACE::lst::kWidth2S
ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kWidth2S
Definition: Common.h:57
ALPAKA_ACCELERATOR_NAMESPACE::lst::CreateQuintuplets::operator()
ALPAKA_FN_ACC void operator()(TAcc const &acc, ModulesConst modules, MiniDoubletsConst mds, SegmentsConst segments, Triplets triplets, TripletsOccupancyConst tripletsOccupancy, Quintuplets quintuplets, QuintupletsOccupancy quintupletsOccupancy, ObjectRangesConst ranges, uint16_t nEligibleT5Modules) const
Definition: Quintuplet.h:2336
ALPAKA_ACCELERATOR_NAMESPACE::lst::computeChiSquared
ALPAKA_FN_ACC ALPAKA_FN_INLINE float computeChiSquared(TAcc const &acc, unsigned int nPoints, float *xs, float *ys, float *delta1, float *delta2, float *slopes, bool *isFlat, float g, float f, float radius)
Definition: Quintuplet.h:1056
ALPAKA_ACCELERATOR_NAMESPACE::lst::T5HasCommonMiniDoublet
ALPAKA_FN_ACC ALPAKA_FN_INLINE bool T5HasCommonMiniDoublet(TripletsConst triplets, SegmentsConst segments, unsigned int innerTripletIndex, unsigned int outerTripletIndex)
Definition: Quintuplet.h:611
ALPAKA_ACCELERATOR_NAMESPACE::lst::kVerticalModuleSlope
ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kVerticalModuleSlope
Definition: Common.h:61
g
The Signals That Services Can Subscribe To This is based on ActivityRegistry and is current per Services can connect to the signals distributed by the ActivityRegistry in order to monitor the activity of the application Each possible callback has some defined which we here list in angle e g
Definition: Activities.doc:4
DiDispStaMuonMonitor_cfi.pt
pt
Definition: DiDispStaMuonMonitor_cfi.py:39
__H2F
#define __H2F
Definition: Common.h:50
__F2H
#define __F2H
Definition: Common.h:49
ALPAKA_ACCELERATOR_NAMESPACE::lst::ptCut
ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float ptCut
Definition: Common.h:52
ALPAKA_ACCELERATOR_NAMESPACE::lst::matchRadiiBBBBE
ALPAKA_FN_ACC ALPAKA_FN_INLINE bool matchRadiiBBBBE(TAcc const &acc, float innerRadius, float bridgeRadius, float outerRadius)
Definition: Quintuplet.h:707
interestingEgammaIsoDetIdsSequence_cff.innerRadius
innerRadius
Definition: interestingEgammaIsoDetIdsSequence_cff.py:15
lst::Quintuplets
QuintupletsSoA::View Quintuplets
Definition: QuintupletsSoA.h:34
ALPAKA_ACCELERATOR_NAMESPACE::lst::matchRadiiBBBEE12378
ALPAKA_FN_ACC ALPAKA_FN_INLINE bool matchRadiiBBBEE12378(TAcc const &acc, float innerRadius, float bridgeRadius, float outerRadius, float bridgeRadiusMin2S, float bridgeRadiusMax2S)
Definition: Quintuplet.h:658
interestingEgammaIsoDetIdsSequence_cff.outerRadius
outerRadius
Definition: interestingEgammaIsoDetIdsSequence_cff.py:14
lst::TripletsConst
TripletsSoA::ConstView TripletsConst
Definition: TripletsSoA.h:31
CosmicsPD_Skims.radius
radius
Definition: CosmicsPD_Skims.py:135
ALPAKA_ACCELERATOR_NAMESPACE::lst::phi_mpi_pi
ALPAKA_FN_HOST_ACC ALPAKA_FN_INLINE float phi_mpi_pi(TAcc const &acc, float x)
Definition: Hit.h:19
ALPAKA_ACCELERATOR_NAMESPACE::lst::runQuintupletDefaultAlgoEEEE
ALPAKA_FN_ACC ALPAKA_FN_INLINE bool runQuintupletDefaultAlgoEEEE(TAcc const &acc, ModulesConst modules, MiniDoubletsConst mds, SegmentsConst segments, uint16_t innerInnerLowerModuleIndex, uint16_t innerOuterLowerModuleIndex, uint16_t outerInnerLowerModuleIndex, uint16_t outerOuterLowerModuleIndex, unsigned int innerSegmentIndex, unsigned int outerSegmentIndex, unsigned int firstMDIndex, unsigned int secondMDIndex, unsigned int thirdMDIndex, unsigned int fourthMDIndex)
Definition: Quintuplet.h:1688
ALPAKA_ACCELERATOR_NAMESPACE::lst::phi
ALPAKA_FN_HOST_ACC ALPAKA_FN_INLINE float phi(TAcc const &acc, float x, float y)
Definition: Hit.h:29
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:23
or
The Signals That Services Can Subscribe To This is based on ActivityRegistry and is current per Services can connect to the signals distributed by the ActivityRegistry in order to monitor the activity of the application Each possible callback has some defined which we here list in angle e< void, edm::EventID const &, edm::Timestamp const & > We also list in braces which AR_WATCH_USING_METHOD_ is used for those or
Definition: Activities.doc:12
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
ALPAKA_ACCELERATOR_NAMESPACE::lst::AddQuintupletRangesToEventExplicit
Definition: Quintuplet.h:2566
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
hcaldqm::ModuleType
ModuleType
Definition: DQModule.h:24
electrons_cff.dzErr
dzErr
Definition: electrons_cff.py:340
f
double f[11][100]
Definition: MuScleFitUtils.cc:78
ALPAKA_ACCELERATOR_NAMESPACE::lst::runDeltaBetaIterationsT5
ALPAKA_FN_ACC ALPAKA_FN_INLINE void runDeltaBetaIterationsT5(TAcc const &acc, float &betaIn, float &betaOut, float betaAv, float &pt_beta, float sdIn_dr, float sdOut_dr, float dr, float lIn)
Definition: Quintuplet.h:1102
ALPAKA_ACCELERATOR_NAMESPACE::lst::t5dnn::kLSTWp2
ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kLSTWp2
Definition: Common.h:67
lst::MiniDoubletsConst
MiniDoubletsSoA::ConstView MiniDoubletsConst
Definition: MiniDoubletsSoA.h:48
ALPAKA_ACCELERATOR_NAMESPACE::lst::runQuintupletDefaultAlgoBBBB
ALPAKA_FN_ACC ALPAKA_FN_INLINE bool runQuintupletDefaultAlgoBBBB(TAcc const &acc, ModulesConst modules, MiniDoubletsConst mds, SegmentsConst segments, uint16_t innerInnerLowerModuleIndex, uint16_t innerOuterLowerModuleIndex, uint16_t outerInnerLowerModuleIndex, uint16_t outerOuterLowerModuleIndex, unsigned int innerSegmentIndex, unsigned int outerSegmentIndex, unsigned int firstMDIndex, unsigned int secondMDIndex, unsigned int thirdMDIndex, unsigned int fourthMDIndex)
Definition: Quintuplet.h:1201
ALPAKA_ACCELERATOR_NAMESPACE::lst::matchRadiiBBEEE
ALPAKA_FN_ACC ALPAKA_FN_INLINE bool matchRadiiBBEEE(TAcc const &acc, float innerRadius, float bridgeRadius, float outerRadius, float bridgeRadiusMin2S, float bridgeRadiusMax2S)
Definition: Quintuplet.h:779
ALCARECOTkAlJpsiMuMu_cff.charge
charge
Definition: ALCARECOTkAlJpsiMuMu_cff.py:47
ALPAKA_ACCELERATOR_NAMESPACE::lst::matchRadiiBBBEE23478
ALPAKA_FN_ACC ALPAKA_FN_INLINE bool matchRadiiBBBEE23478(TAcc const &acc, float innerRadius, float bridgeRadius, float outerRadius, float bridgeRadiusMin2S, float bridgeRadiusMax2S)
Definition: Quintuplet.h:731
correctionTermsCaloMet_cff.C
C
Definition: correctionTermsCaloMet_cff.py:34
lst::ModulesConst
ModulesSoA::ConstView ModulesConst
Definition: ModulesSoA.h:47
ALPAKA_ACCELERATOR_NAMESPACE::lst::passT5RZConstraint
ALPAKA_FN_ACC ALPAKA_FN_INLINE bool passT5RZConstraint(TAcc const &acc, ModulesConst modules, MiniDoubletsConst mds, unsigned int firstMDIndex, unsigned int secondMDIndex, unsigned int thirdMDIndex, unsigned int fourthMDIndex, unsigned int fifthMDIndex, uint16_t lowerModuleIndex1, uint16_t lowerModuleIndex2, uint16_t lowerModuleIndex3, uint16_t lowerModuleIndex4, uint16_t lowerModuleIndex5, float &rzChiSquared, float inner_pt, float innerRadius, float g, float f, bool &tightCutFlag)
Definition: Quintuplet.h:179
ALPAKA_ACCELERATOR_NAMESPACE::lst::matchRadiiBBBBB
ALPAKA_FN_ACC ALPAKA_FN_INLINE bool matchRadiiBBBBB(TAcc const &acc, float innerRadius, float bridgeRadius, float outerRadius)
Definition: Quintuplet.h:683
tas::module_layers
const std::vector< int > & module_layers()
Definition: LSTEff.cc:3388
lst::QuintupletsOccupancy
QuintupletsOccupancySoA::View QuintupletsOccupancy
Definition: QuintupletsSoA.h:42
diffTwoXMLs.ranges
string ranges
Definition: diffTwoXMLs.py:79
ALPAKA_ACCELERATOR_NAMESPACE::lst::computeErrorInRadius
ALPAKA_FN_ACC ALPAKA_FN_INLINE void computeErrorInRadius(TAcc const &acc, float *x1Vec, float *y1Vec, float *x2Vec, float *y2Vec, float *x3Vec, float *y3Vec, float &minimumRadius, float &maximumRadius)
Definition: Quintuplet.h:626
ALPAKA_ACCELERATOR_NAMESPACE::lst::kMulsInGeV
ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kMulsInGeV
Definition: Common.h:41
ALPAKA_ACCELERATOR_NAMESPACE::lst::kPixelPSZpitch
ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kPixelPSZpitch
Definition: Common.h:54
ALPAKA_ACCELERATOR_NAMESPACE::lst::matchRadiiEEEEE
ALPAKA_FN_ACC ALPAKA_FN_INLINE bool matchRadiiEEEEE(TAcc const &acc, float innerRadius, float bridgeRadius, float outerRadius, float innerRadiusMin2S, float innerRadiusMax2S, float bridgeRadiusMin2S, float bridgeRadiusMax2S)
Definition: Quintuplet.h:840
tas::module_subdets
const std::vector< int > & module_subdets()
Definition: LSTEff.cc:3288
ALPAKA_ACCELERATOR_NAMESPACE::lst::kStripPSZpitch
ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kStripPSZpitch
Definition: Common.h:55
ALPAKA_ACCELERATOR_NAMESPACE::lst::kPt_betaMax
ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kPt_betaMax
Definition: Common.h:59
ALPAKA_ACCELERATOR_NAMESPACE::lst::kPi
ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kPi
Definition: Common.h:39
ALPAKA_ACCELERATOR_NAMESPACE::lst::runQuintupletDefaultAlgoBBEE
ALPAKA_FN_ACC ALPAKA_FN_INLINE bool runQuintupletDefaultAlgoBBEE(TAcc const &acc, ModulesConst modules, MiniDoubletsConst mds, SegmentsConst segments, uint16_t innerInnerLowerModuleIndex, uint16_t innerOuterLowerModuleIndex, uint16_t outerInnerLowerModuleIndex, uint16_t outerOuterLowerModuleIndex, unsigned int innerSegmentIndex, unsigned int outerSegmentIndex, unsigned int firstMDIndex, unsigned int secondMDIndex, unsigned int thirdMDIndex, unsigned int fourthMDIndex)
Definition: Quintuplet.h:1443
ALPAKA_ACCELERATOR_NAMESPACE::lst::kStrip2SZpitch
ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kStrip2SZpitch
Definition: Common.h:56
lst::Triplets
TripletsSoA::View Triplets
Definition: TripletsSoA.h:30
testProducerWithPsetDescEmpty_cfi.y1
y1
Definition: testProducerWithPsetDescEmpty_cfi.py:29
ALPAKA_ACCELERATOR_NAMESPACE::lst::CreateEligibleModulesListForQuintuplets::operator()
ALPAKA_FN_ACC void operator()(TAcc const &acc, ModulesConst modules, TripletsOccupancyConst tripletsOccupancy, ObjectRanges ranges) const
Definition: Quintuplet.h:2464
ALPAKA_ACCELERATOR_NAMESPACE::lst::AddQuintupletRangesToEventExplicit::operator()
ALPAKA_FN_ACC void operator()(TAcc const &acc, ModulesConst modules, QuintupletsOccupancyConst quintupletsOccupancy, ObjectRanges ranges) const
Definition: Quintuplet.h:2568
lst::TripletsOccupancyConst
TripletsOccupancySoA::ConstView TripletsOccupancyConst
Definition: TripletsSoA.h:39
lst::ObjectRangesConst
ObjectRangesSoA::ConstView ObjectRangesConst
Definition: ObjectRangesSoA.h:32
a
double a
Definition: hdecay.h:121
ALPAKA_ACCELERATOR_NAMESPACE::lst::kSinAlphaMax
ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kSinAlphaMax
Definition: Common.h:51
lst::SegmentsConst
SegmentsSoA::ConstView SegmentsConst
Definition: SegmentsSoA.h:49
ppsModifySingularModes_cfi.z1
z1
Definition: ppsModifySingularModes_cfi.py:8
x
float x
Definition: beamSpotDipStandalone.cc:55
lst::QuintupletsOccupancyConst
QuintupletsOccupancySoA::ConstView QuintupletsOccupancyConst
Definition: QuintupletsSoA.h:43
ALPAKA_ACCELERATOR_NAMESPACE::lst::matchRadiiBEEEE
ALPAKA_FN_ACC ALPAKA_FN_INLINE bool matchRadiiBEEEE(TAcc const &acc, float innerRadius, float bridgeRadius, float outerRadius, float innerRadiusMin2S, float innerRadiusMax2S, float bridgeRadiusMin2S, float bridgeRadiusMax2S)
Definition: Quintuplet.h:809
A
Definition: APVGainStruct.h:7
ALPAKA_ACCELERATOR_NAMESPACE::lst::addQuintupletToMemory
ALPAKA_FN_ACC ALPAKA_FN_INLINE void addQuintupletToMemory(TripletsConst triplets, Quintuplets quintuplets, unsigned int innerTripletIndex, unsigned int outerTripletIndex, uint16_t lowerModule1, uint16_t lowerModule2, uint16_t lowerModule3, uint16_t lowerModule4, uint16_t lowerModule5, float innerRadius, float bridgeRadius, float outerRadius, float regressionG, float regressionF, float regressionRadius, float rzChiSquared, float rPhiChiSquared, float nonAnchorChiSquared, float pt, float eta, float phi, float scores, uint8_t layer, unsigned int quintupletIndex, bool tightCutFlag)
Definition: Quintuplet.h:28
ALPAKA_ACCELERATOR_NAMESPACE::lst::kDeltaZLum
ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kDeltaZLum
Definition: Common.h:53
ALPAKA_ACCELERATOR_NAMESPACE::lst::deltaPhi
ALPAKA_FN_HOST_ACC ALPAKA_FN_INLINE float deltaPhi(TAcc const &acc, float x1, float y1, float x2, float y2)
Definition: Hit.h:34
ALPAKA_ACCELERATOR_NAMESPACE::caPixelDoublets::ALPAKA_ASSERT_ACC
ALPAKA_ASSERT_ACC(offsets)
fireworks::subdets
static const std::string subdets[7]
Definition: TrackUtils.cc:60
cms::cudacompat::atomicAdd
T1 atomicAdd(T1 *a, T2 b)
Definition: cudaCompat.h:61
AlCaHLTBitMon_ParallelJobs.p
def p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
ALPAKA_ACCELERATOR_NAMESPACE::lst::computeRadiusUsingRegression
ALPAKA_FN_ACC ALPAKA_FN_INLINE float computeRadiusUsingRegression(TAcc const &acc, unsigned int nPoints, float *xs, float *ys, float *delta1, float *delta2, float *slopes, bool *isFlat, float &g, float &f, float *sigmas2, float &chiSquared)
Definition: Quintuplet.h:958
ALPAKA_ACCELERATOR_NAMESPACE::lst::matchRadiiBBBEE34578
ALPAKA_FN_ACC ALPAKA_FN_INLINE bool matchRadiiBBBEE34578(TAcc const &acc, float innerRadius, float bridgeRadius, float outerRadius, float bridgeRadiusMin2S, float bridgeRadiusMax2S)
Definition: Quintuplet.h:755
ALPAKA_ACCELERATOR_NAMESPACE::lst::CreateEligibleModulesListForQuintuplets
Definition: Quintuplet.h:2462
ALPAKA_ACCELERATOR_NAMESPACE::lst::eta
ALPAKA_FN_HOST_ACC ALPAKA_FN_INLINE float eta(TAcc const &acc, float x, float y, float z)
Definition: Hit.h:11
Generated for CMSSW Reference Manual by doxygen 1.8.14