1 #ifndef RecoTracker_LSTCore_src_alpaka_PixelTriplet_h 2 #define RecoTracker_LSTCore_src_alpaka_PixelTriplet_h 15 template <
typename TAcc>
22 uint16_t pixelModuleIndex,
23 uint16_t outerInnerLowerModuleIndex,
24 uint16_t outerOuterLowerModuleIndex,
25 unsigned int innerSegmentIndex,
26 unsigned int outerSegmentIndex,
27 unsigned int firstMDIndex,
28 unsigned int secondMDIndex,
29 unsigned int thirdMDIndex,
30 unsigned int fourthMDIndex);
31 template <
typename TAcc>
38 uint16_t pixelModuleIndex,
39 uint16_t outerInnerLowerModuleIndex,
40 uint16_t outerOuterLowerModuleIndex,
41 unsigned int innerSegmentIndex,
42 unsigned int outerSegmentIndex,
43 unsigned int firstMDIndex,
44 unsigned int secondMDIndex,
45 unsigned int thirdMDIndex,
46 unsigned int fourthMDIndex);
52 unsigned int pixelSegmentIndex,
53 unsigned int tripletIndex,
59 float rPhiChiSquaredInwards,
61 unsigned int pixelTripletIndex,
68 pixelTriplets.pixelSegmentIndices()[pixelTripletIndex] = pixelSegmentIndex;
69 pixelTriplets.tripletIndices()[pixelTripletIndex] = tripletIndex;
70 pixelTriplets.pixelRadius()[pixelTripletIndex] =
__F2H(pixelRadius);
71 pixelTriplets.tripletRadius()[pixelTripletIndex] =
__F2H(tripletRadius);
72 pixelTriplets.pt()[pixelTripletIndex] =
__F2H(
pt);
73 pixelTriplets.eta()[pixelTripletIndex] =
__F2H(
eta);
74 pixelTriplets.phi()[pixelTripletIndex] =
__F2H(
phi);
75 pixelTriplets.eta_pix()[pixelTripletIndex] =
__F2H(eta_pix);
76 pixelTriplets.phi_pix()[pixelTripletIndex] =
__F2H(phi_pix);
77 pixelTriplets.isDup()[pixelTripletIndex] =
false;
78 pixelTriplets.score()[pixelTripletIndex] =
__F2H(
score);
80 pixelTriplets.centerX()[pixelTripletIndex] =
__F2H(centerX);
81 pixelTriplets.centerY()[pixelTripletIndex] =
__F2H(centerY);
82 pixelTriplets.logicalLayers()[pixelTripletIndex][0] = 0;
83 pixelTriplets.logicalLayers()[pixelTripletIndex][1] = 0;
84 pixelTriplets.logicalLayers()[pixelTripletIndex][2] = triplets.logicalLayers()[tripletIndex][0];
85 pixelTriplets.logicalLayers()[pixelTripletIndex][3] = triplets.logicalLayers()[tripletIndex][1];
86 pixelTriplets.logicalLayers()[pixelTripletIndex][4] = triplets.logicalLayers()[tripletIndex][2];
88 pixelTriplets.lowerModuleIndices()[pixelTripletIndex][0] = segments.innerLowerModuleIndices()[pixelSegmentIndex];
89 pixelTriplets.lowerModuleIndices()[pixelTripletIndex][1] = segments.outerLowerModuleIndices()[pixelSegmentIndex];
90 pixelTriplets.lowerModuleIndices()[pixelTripletIndex][2] = triplets.lowerModuleIndices()[tripletIndex][0];
91 pixelTriplets.lowerModuleIndices()[pixelTripletIndex][3] = triplets.lowerModuleIndices()[tripletIndex][1];
92 pixelTriplets.lowerModuleIndices()[pixelTripletIndex][4] = triplets.lowerModuleIndices()[tripletIndex][2];
94 unsigned int pixelInnerMD = segments.mdIndices()[pixelSegmentIndex][0];
95 unsigned int pixelOuterMD = segments.mdIndices()[pixelSegmentIndex][1];
97 pixelTriplets.hitIndices()[pixelTripletIndex][0] = mds.anchorHitIndices()[pixelInnerMD];
98 pixelTriplets.hitIndices()[pixelTripletIndex][1] = mds.outerHitIndices()[pixelInnerMD];
99 pixelTriplets.hitIndices()[pixelTripletIndex][2] = mds.anchorHitIndices()[pixelOuterMD];
100 pixelTriplets.hitIndices()[pixelTripletIndex][3] = mds.outerHitIndices()[pixelOuterMD];
102 pixelTriplets.hitIndices()[pixelTripletIndex][4] = triplets.hitIndices()[tripletIndex][0];
103 pixelTriplets.hitIndices()[pixelTripletIndex][5] = triplets.hitIndices()[tripletIndex][1];
104 pixelTriplets.hitIndices()[pixelTripletIndex][6] = triplets.hitIndices()[tripletIndex][2];
105 pixelTriplets.hitIndices()[pixelTripletIndex][7] = triplets.hitIndices()[tripletIndex][3];
106 pixelTriplets.hitIndices()[pixelTripletIndex][8] = triplets.hitIndices()[tripletIndex][4];
107 pixelTriplets.hitIndices()[pixelTripletIndex][9] = triplets.hitIndices()[tripletIndex][5];
108 pixelTriplets.rPhiChiSquared()[pixelTripletIndex] = rPhiChiSquared;
109 pixelTriplets.rPhiChiSquaredInwards()[pixelTripletIndex] = rPhiChiSquaredInwards;
110 pixelTriplets.rzChiSquared()[pixelTripletIndex] = rzChiSquared;
113 template <
typename TAcc>
120 uint16_t pixelLowerModuleIndex,
121 uint16_t outerInnerLowerModuleIndex,
122 uint16_t outerOuterLowerModuleIndex,
123 unsigned int innerSegmentIndex,
124 unsigned int outerSegmentIndex) {
125 short outerInnerLowerModuleSubdet =
modules.subdets()[outerInnerLowerModuleIndex];
126 short outerOuterLowerModuleSubdet =
modules.subdets()[outerOuterLowerModuleIndex];
128 unsigned int firstMDIndex = segments.mdIndices()[innerSegmentIndex][0];
129 unsigned int secondMDIndex = segments.mdIndices()[innerSegmentIndex][1];
131 unsigned int thirdMDIndex = segments.mdIndices()[outerSegmentIndex][0];
132 unsigned int fourthMDIndex = segments.mdIndices()[outerSegmentIndex][1];
134 if (outerInnerLowerModuleSubdet == Barrel and
135 (outerOuterLowerModuleSubdet == Barrel
or outerOuterLowerModuleSubdet ==
Endcap)) {
142 pixelLowerModuleIndex,
143 outerInnerLowerModuleIndex,
144 outerOuterLowerModuleIndex,
151 }
else if (outerInnerLowerModuleSubdet ==
Endcap and outerOuterLowerModuleSubdet ==
Endcap) {
158 pixelLowerModuleIndex,
159 outerInnerLowerModuleIndex,
160 outerOuterLowerModuleIndex,
172 uint16_t lowerModuleIndex1,
173 uint16_t lowerModuleIndex2,
174 uint16_t lowerModuleIndex3,
175 float rzChiSquared) {
186 if (layer1 == 8 and layer2 == 9 and layer3 == 10) {
187 return rzChiSquared < 13.6067f;
188 }
else if (layer1 == 8 and layer2 == 9 and layer3 == 15) {
189 return rzChiSquared < 5.5953f;
190 }
else if (layer1 == 7 and layer2 == 8 and layer3 == 9) {
191 return rzChiSquared < 3.9263f;
200 else if (layer1 == 1 and layer2 == 2 and layer3 == 3) {
201 return rzChiSquared < 9.4377f;
202 }
else if (layer1 == 1 and layer2 == 2 and layer3 == 7) {
203 return rzChiSquared < 9.9975f;
204 }
else if (layer1 == 1 and layer2 == 7 and layer3 == 8) {
205 return rzChiSquared < 8.6369f;
206 }
else if (layer1 == 2 and layer2 == 3 and layer3 == 7) {
207 return rzChiSquared < 37.945f;
208 }
else if (layer1 == 2 and layer2 == 3 and layer3 == 12) {
209 return rzChiSquared < 43.0167f;
210 }
else if (layer1 == 2 and layer2 == 3 and layer3 == 4) {
211 return rzChiSquared < 8.6923f;
212 }
else if (layer1 == 2 and layer2 == 7 and layer3 == 8) {
213 return rzChiSquared < 11.9672f;
214 }
else if (layer1 == 2 and layer2 == 7 and layer3 == 13) {
215 return rzChiSquared < 16.2133f;
222 template <
typename TAcc>
224 unsigned int nPoints,
237 float chiSquared = 0.f;
238 float absArctanSlope, angleM, xPrime, yPrime, sigma2;
239 for (
size_t i = 0;
i < nPoints;
i++) {
242 if (xs[
i] > 0 and ys[
i] > 0) {
243 angleM =
kPi / 2.f - absArctanSlope;
244 }
else if (xs[
i] < 0 and ys[
i] > 0) {
245 angleM = absArctanSlope +
kPi / 2.f;
246 }
else if (xs[
i] < 0 and ys[
i] < 0) {
247 angleM = -(absArctanSlope +
kPi / 2.f);
248 }
else if (xs[
i] > 0 and ys[
i] < 0) {
249 angleM = -(
kPi / 2.f - absArctanSlope);
261 sigma2 = 4 * ((xPrime * delta1[
i]) * (xPrime * delta1[
i]) + (yPrime * delta2[
i]) * (yPrime * delta2[
i]));
262 chiSquared += (xs[
i] * xs[
i] + ys[
i] * ys[
i] - 2 *
g * xs[
i] - 2 *
f * ys[
i] +
c) *
263 (xs[
i] * xs[
i] + ys[
i] * ys[
i] - 2 *
g * xs[
i] - 2 *
f * ys[
i] +
c) / sigma2;
269 template <
typename TAcc>
272 uint16_t* lowerModuleIndices,
278 float delta1[3]{}, delta2[3]{}, slopes[3]{};
280 float chiSquared = 0;
283 for (
size_t i = 0;
i < 3;
i++) {
285 short moduleSubdet =
modules.subdets()[lowerModuleIndices[
i]];
286 short moduleSide =
modules.sides()[lowerModuleIndices[
i]];
287 float drdz =
modules.drdzs()[lowerModuleIndices[
i]];
288 slopes[
i] =
modules.dxdys()[lowerModuleIndices[
i]];
290 if (moduleSubdet == Barrel and moduleType ==
PS and moduleSide ==
Center) {
297 else if (moduleSubdet == Barrel and moduleType ==
TwoS) {
304 else if (moduleSubdet == Barrel and moduleType ==
PS and moduleSide !=
Center) {
310 else if (moduleSubdet ==
Endcap and moduleType ==
PS) {
321 else if (moduleSubdet ==
Endcap and moduleType ==
TwoS) {
323 delta2[
i] = 500 * inv1;
328 printf(
"ERROR!!!!! I SHOULDN'T BE HERE!!!! subdet = %d, type = %d, side = %d\n",
335 chiSquared =
computeChiSquaredpT3(acc, 3, xs, ys, delta1, delta2, slopes, isFlat,
g,
f,
radius);
341 float g,
float f,
float r,
float* xPix,
float* yPix) {
342 float residual = (xPix[0] -
g) * (xPix[0] -
g) + (yPix[0] -
f) * (yPix[0] -
f) - r * r;
343 float chiSquared = residual * residual;
344 residual = (xPix[1] -
g) * (xPix[1] -
g) + (yPix[1] -
f) * (yPix[1] -
f) - r * r;
345 chiSquared += residual * residual;
353 uint16_t lowerModuleIndex1,
354 uint16_t lowerModuleIndex2,
355 uint16_t lowerModuleIndex3,
367 if (layer1 == 8 and layer2 == 9 and layer3 == 10) {
368 return chiSquared < 7.003f;
369 }
else if (layer1 == 8 and layer2 == 9 and layer3 == 15) {
370 return chiSquared < 0.5f;
371 }
else if (layer1 == 7 and layer2 == 8 and layer3 == 9) {
372 return chiSquared < 8.046f;
373 }
else if (layer1 == 7 and layer2 == 8 and layer3 == 14) {
374 return chiSquared < 0.575f;
375 }
else if (layer1 == 1 and layer2 == 2 and layer3 == 7) {
376 return chiSquared < 5.304f;
377 }
else if (layer1 == 1 and layer2 == 2 and layer3 == 3) {
378 return chiSquared < 10.6211f;
379 }
else if (layer1 == 1 and layer2 == 7 and layer3 == 8) {
380 return chiSquared < 4.617f;
381 }
else if (layer1 == 2 and layer2 == 7 and layer3 == 8) {
382 return chiSquared < 8.046f;
383 }
else if (layer1 == 2 and layer2 == 7 and layer3 == 13) {
384 return chiSquared < 0.435f;
385 }
else if (layer1 == 2 and layer2 == 3 and layer3 == 7) {
386 return chiSquared < 9.244f;
387 }
else if (layer1 == 2 and layer2 == 3 and layer3 == 12) {
388 return chiSquared < 0.287f;
389 }
else if (layer1 == 2 and layer2 == 3 and layer3 == 4) {
390 return chiSquared < 18.509f;
397 uint16_t lowerModuleIndex1,
398 uint16_t lowerModuleIndex2,
399 uint16_t lowerModuleIndex3,
411 if (layer1 == 7 and layer2 == 8 and layer3 == 9)
413 return chiSquared < 22016.8055f;
414 }
else if (layer1 == 7 and layer2 == 8 and layer3 == 14)
416 return chiSquared < 935179.56807f;
417 }
else if (layer1 == 8 and layer2 == 9 and layer3 == 10)
419 return chiSquared < 29064.12959f;
420 }
else if (layer1 == 8 and layer2 == 9 and layer3 == 15)
422 return chiSquared < 935179.5681f;
423 }
else if (layer1 == 1 and layer2 == 2 and layer3 == 3)
425 return chiSquared < 1370.0113195101474f;
426 }
else if (layer1 == 1 and layer2 == 2 and layer3 == 7)
428 return chiSquared < 5492.110048314815f;
429 }
else if (layer1 == 2 and layer2 == 3 and layer3 == 4)
431 return chiSquared < 4160.410806470067f;
432 }
else if (layer1 == 1 and layer2 == 7 and layer3 == 8)
434 return chiSquared < 29064.129591225726f;
435 }
else if (layer1 == 2 and layer2 == 3 and layer3 == 7)
437 return chiSquared < 12634.215376250893f;
438 }
else if (layer1 == 2 and layer2 == 3 and layer3 == 12)
440 return chiSquared < 353821.69361145404f;
441 }
else if (layer1 == 2 and layer2 == 7 and layer3 == 8)
443 return chiSquared < 33393.26076341235f;
444 }
else if (layer1 == 2 and layer2 == 7 and layer3 == 13)
446 return chiSquared < 935179.5680742573f;
456 return ((firstMin <= secondMin) && (secondMin < firstMax)) || ((secondMin < firstMin) && (firstMin < secondMax));
460 template <
typename TAcc>
463 float pixelRadiusError,
464 float tripletRadius) {
465 float tripletInvRadiusErrorBound = 0.15624f;
466 float pixelInvRadiusErrorBound = 0.17235f;
469 pixelInvRadiusErrorBound = 0.6375f;
470 tripletInvRadiusErrorBound = 0.6588f;
473 float tripletRadiusInvMax = (1 + tripletInvRadiusErrorBound) / tripletRadius;
474 float tripletRadiusInvMin =
alpaka::math::max(acc, (1 - tripletInvRadiusErrorBound) / tripletRadius, 0.0
f);
476 float pixelRadiusInvMax =
477 alpaka::math::max(acc, (1 + pixelInvRadiusErrorBound) / pixelRadius, 1.
f / (pixelRadius - pixelRadiusError));
478 float pixelRadiusInvMin =
479 alpaka::math::min(acc, (1 - pixelInvRadiusErrorBound) / pixelRadius, 1.
f / (pixelRadius + pixelRadiusError));
481 return checkIntervalOverlappT3(tripletRadiusInvMin, tripletRadiusInvMax, pixelRadiusInvMin, pixelRadiusInvMax);
484 template <
typename TAcc>
487 float pixelRadiusError,
488 float tripletRadius) {
489 float tripletInvRadiusErrorBound = 0.45972f;
490 float pixelInvRadiusErrorBound = 0.19644f;
493 pixelInvRadiusErrorBound = 0.6805f;
494 tripletInvRadiusErrorBound = 0.8557f;
497 float tripletRadiusInvMax = (1 + tripletInvRadiusErrorBound) / tripletRadius;
498 float tripletRadiusInvMin =
alpaka::math::max(acc, (1 - tripletInvRadiusErrorBound) / tripletRadius, 0.0
f);
500 float pixelRadiusInvMax =
501 alpaka::math::max(acc, (1 + pixelInvRadiusErrorBound) / pixelRadius, 1.
f / (pixelRadius - pixelRadiusError));
502 float pixelRadiusInvMin =
503 alpaka::math::min(acc, (1 - pixelInvRadiusErrorBound) / pixelRadius, 1.
f / (pixelRadius + pixelRadiusError));
505 return checkIntervalOverlappT3(tripletRadiusInvMin, tripletRadiusInvMax, pixelRadiusInvMin, pixelRadiusInvMax);
508 template <
typename TAcc>
511 float pixelRadiusError,
512 float tripletRadius) {
513 float tripletInvRadiusErrorBound = 1.59294f;
514 float pixelInvRadiusErrorBound = 0.255181f;
518 pixelInvRadiusErrorBound = 2.2091f;
519 tripletInvRadiusErrorBound = 2.3548f;
522 float tripletRadiusInvMax = (1 + tripletInvRadiusErrorBound) / tripletRadius;
523 float tripletRadiusInvMin =
alpaka::math::max(acc, (1 - tripletInvRadiusErrorBound) / tripletRadius, 0.0
f);
525 float pixelRadiusInvMax =
526 alpaka::math::max(acc, (1 + pixelInvRadiusErrorBound) / pixelRadius, 1.
f / (pixelRadius - pixelRadiusError));
527 float pixelRadiusInvMin =
528 alpaka::math::min(acc, (1 - pixelInvRadiusErrorBound) / pixelRadius, 1.
f / (pixelRadius + pixelRadiusError));
531 return checkIntervalOverlappT3(tripletRadiusInvMin, tripletRadiusInvMax, pixelRadiusInvMin, pixelRadiusInvMax);
534 template <
typename TAcc>
537 float pixelRadiusError,
538 float tripletRadius) {
539 float tripletInvRadiusErrorBound = 1.7006f;
540 float pixelInvRadiusErrorBound = 0.26367f;
544 pixelInvRadiusErrorBound = 2.286f;
545 tripletInvRadiusErrorBound = 2.436f;
548 float tripletRadiusInvMax = (1 + tripletInvRadiusErrorBound) / tripletRadius;
549 float tripletRadiusInvMin =
alpaka::math::max(acc, (1 - tripletInvRadiusErrorBound) / tripletRadius, 0.0
f);
551 float pixelRadiusInvMax =
552 alpaka::math::max(acc, (1 + pixelInvRadiusErrorBound) / pixelRadius, 1.
f / (pixelRadius - pixelRadiusError));
553 float pixelRadiusInvMin =
554 alpaka::math::min(acc, (1 - pixelInvRadiusErrorBound) / pixelRadius, 1.
f / (pixelRadius + pixelRadiusError));
557 return checkIntervalOverlappT3(tripletRadiusInvMin, tripletRadiusInvMax, pixelRadiusInvMin, pixelRadiusInvMax);
560 template <
typename TAcc>
564 float pixelRadiusError,
566 int16_t lowerModuleIndex,
567 uint16_t middleModuleIndex,
568 uint16_t upperModuleIndex) {
580 template <
typename TAcc>
583 const uint16_t* lowerModuleIndices,
592 float pixelSegmentPt,
593 float pixelSegmentPx,
594 float pixelSegmentPy,
595 float pixelSegmentPz,
596 int pixelSegmentCharge) {
601 float Px = pixelSegmentPx, Py = pixelSegmentPy, Pz = pixelSegmentPz;
602 int charge = pixelSegmentCharge;
603 float x1 = xPix[1] / 100;
604 float y1 = yPix[1] / 100;
605 float z1 = zPix[1] / 100;
606 float r1 = rtPix[1] / 100;
610 for (
size_t i = 0;
i < Params_T3::kLayers;
i++) {
611 float zsi = zs[
i] / 100;
612 float rtsi = rts[
i] / 100;
613 uint16_t lowerModuleIndex = lowerModuleIndices[
i];
614 const int moduleType =
modules.moduleType()[lowerModuleIndex];
615 const int moduleSide =
modules.sides()[lowerModuleIndex];
616 const int moduleSubdet =
modules.subdets()[lowerModuleIndex];
623 if (moduleSubdet ==
Endcap) {
624 float s = (zsi -
z1) *
p / Pz;
630 if (moduleSubdet == Barrel) {
631 float paraA = r1 * r1 + 2 * (Px * Px + Py * Py) / (
a *
a) + 2 * (
y1 * Px - x1 * Py) /
a - rtsi * rtsi;
632 float paraB = 2 * (x1 * Px +
y1 * Py) /
a;
633 float paraC = 2 * (
y1 * Px - x1 * Py) /
a + 2 * (Px * Px + Py * Py) / (
a *
a);
634 float A = paraB * paraB + paraC * paraC;
635 float B = 2 * paraA * paraB;
636 float C = paraA * paraA - paraC * paraC;
639 float solz1 = alpaka::math::asin(acc, sol1) / rou * Pz /
p +
z1;
640 float solz2 = alpaka::math::asin(acc, sol2) / rou * Pz /
p +
z1;
646 residual = moduleSubdet == Barrel ? diffz : diffr;
649 if (moduleType == 0) {
657 if (moduleType == 0 and moduleSubdet == Barrel and moduleSide !=
Center) {
658 float drdz =
modules.drdzs()[lowerModuleIndex];
659 error2 /= (1 + drdz * drdz);
661 RMSE += (residual * residual) / error2;
669 template <
typename TAcc>
677 unsigned int pixelSegmentIndex,
678 unsigned int tripletIndex,
680 float& tripletRadius,
684 float& rPhiChiSquared,
685 float& rPhiChiSquaredInwards,
686 bool runChiSquaredCuts =
true) {
688 uint16_t pixelModuleIndex = segments.innerLowerModuleIndices()[pixelSegmentIndex];
690 uint16_t lowerModuleIndex = triplets.lowerModuleIndices()[tripletIndex][0];
691 uint16_t middleModuleIndex = triplets.lowerModuleIndices()[tripletIndex][1];
692 uint16_t upperModuleIndex = triplets.lowerModuleIndices()[tripletIndex][2];
706 triplets.segmentIndices()[tripletIndex][0]))
720 triplets.segmentIndices()[tripletIndex][1]))
725 unsigned int pixelSegmentArrayIndex = pixelSegmentIndex -
ranges.segmentModuleIndices()[pixelModuleIndex];
726 float pixelSegmentPt = segmentsPixel.ptIn()[pixelSegmentArrayIndex];
727 float pixelSegmentPtError = segmentsPixel.ptErr()[pixelSegmentArrayIndex];
728 float pixelSegmentPx = segmentsPixel.px()[pixelSegmentArrayIndex];
729 float pixelSegmentPy = segmentsPixel.py()[pixelSegmentArrayIndex];
730 float pixelSegmentPz = segmentsPixel.pz()[pixelSegmentArrayIndex];
731 int pixelSegmentCharge = segmentsPixel.charge()[pixelSegmentArrayIndex];
733 float pixelG = segmentsPixel.circleCenterX()[pixelSegmentArrayIndex];
734 float pixelF = segmentsPixel.circleCenterY()[pixelSegmentArrayIndex];
735 float pixelRadiusPCA = segmentsPixel.circleRadius()[pixelSegmentArrayIndex];
737 unsigned int pixelInnerMDIndex = segments.mdIndices()[pixelSegmentIndex][0];
738 unsigned int pixelOuterMDIndex = segments.mdIndices()[pixelSegmentIndex][1];
740 pixelRadius = pixelSegmentPt *
kR1GeVf;
741 float pixelRadiusError = pixelSegmentPtError *
kR1GeVf;
742 unsigned int tripletInnerSegmentIndex = triplets.segmentIndices()[tripletIndex][0];
743 unsigned int tripletOuterSegmentIndex = triplets.segmentIndices()[tripletIndex][1];
745 unsigned int firstMDIndex = segments.mdIndices()[tripletInnerSegmentIndex][0];
746 unsigned int secondMDIndex = segments.mdIndices()[tripletInnerSegmentIndex][1];
747 unsigned int thirdMDIndex = segments.mdIndices()[tripletOuterSegmentIndex][1];
749 float xs[Params_T3::kLayers] = {
750 mds.anchorX()[firstMDIndex], mds.anchorX()[secondMDIndex], mds.anchorX()[thirdMDIndex]};
751 float ys[Params_T3::kLayers] = {
752 mds.anchorY()[firstMDIndex], mds.anchorY()[secondMDIndex], mds.anchorY()[thirdMDIndex]};
755 tripletRadius = triplets.radius()[tripletIndex];
756 g = triplets.centerX()[tripletIndex];
757 f = triplets.centerY()[tripletIndex];
769 uint16_t lowerModuleIndices[Params_T3::kLayers] = {lowerModuleIndex, middleModuleIndex, upperModuleIndex};
771 if (runChiSquaredCuts and pixelSegmentPt < 5.0
f) {
772 float rts[Params_T3::kLayers] = {
773 mds.anchorRt()[firstMDIndex], mds.anchorRt()[secondMDIndex], mds.anchorRt()[thirdMDIndex]};
774 float zs[Params_T3::kLayers] = {
775 mds.anchorZ()[firstMDIndex], mds.anchorZ()[secondMDIndex], mds.anchorZ()[thirdMDIndex]};
776 float rtPix[Params_pLS::kLayers] = {mds.anchorRt()[pixelInnerMDIndex], mds.anchorRt()[pixelOuterMDIndex]};
777 float xPix[Params_pLS::kLayers] = {mds.anchorX()[pixelInnerMDIndex], mds.anchorX()[pixelOuterMDIndex]};
778 float yPix[Params_pLS::kLayers] = {mds.anchorY()[pixelInnerMDIndex], mds.anchorY()[pixelOuterMDIndex]};
779 float zPix[Params_pLS::kLayers] = {mds.anchorZ()[pixelInnerMDIndex], mds.anchorZ()[pixelOuterMDIndex]};
805 if (runChiSquaredCuts and pixelSegmentPt < 5.0
f) {
810 float xPix[Params_pLS::kLayers] = {mds.anchorX()[pixelInnerMDIndex], mds.anchorX()[pixelOuterMDIndex]};
811 float yPix[Params_pLS::kLayers] = {mds.anchorY()[pixelInnerMDIndex], mds.anchorY()[pixelOuterMDIndex]};
814 if (runChiSquaredCuts and pixelSegmentPt < 5.0
f) {
816 modules, lowerModuleIndex, middleModuleIndex, upperModuleIndex, rPhiChiSquaredInwards))
825 template <
typename TAcc>
836 unsigned int* connectedPixelSize,
837 unsigned int* connectedPixelIndex,
838 unsigned int nPixelSegments)
const {
839 auto const globalBlockIdx = alpaka::getIdx<alpaka::Grid, alpaka::Blocks>(acc);
840 auto const globalThreadIdx = alpaka::getIdx<alpaka::Grid, alpaka::Threads>(acc);
841 auto const gridBlockExtent = alpaka::getWorkDiv<alpaka::Grid, alpaka::Blocks>(acc);
842 auto const gridThreadExtent = alpaka::getWorkDiv<alpaka::Grid, alpaka::Threads>(acc);
844 for (
unsigned int i_pLS = globalThreadIdx[1]; i_pLS < nPixelSegments; i_pLS += gridThreadExtent[1]) {
845 auto iLSModule_max = connectedPixelIndex[i_pLS] + connectedPixelSize[i_pLS];
847 for (
unsigned int iLSModule = connectedPixelIndex[i_pLS] + globalBlockIdx[0]; iLSModule < iLSModule_max;
848 iLSModule += gridBlockExtent[0]) {
849 uint16_t tripletLowerModuleIndex =
850 modulesPixel.connectedPixels()
853 if (tripletLowerModuleIndex >=
modules.nLowerModules()) {
854 printf(
"tripletLowerModuleIndex %d >= modules.nLowerModules %d \n",
855 tripletLowerModuleIndex,
861 if (
modules.moduleType()[tripletLowerModuleIndex] ==
TwoS)
864 uint16_t pixelModuleIndex =
modules.nLowerModules();
865 unsigned int nOuterTriplets = tripletsOccupancy.nTriplets()[tripletLowerModuleIndex];
866 if (nOuterTriplets == 0)
869 unsigned int pixelSegmentIndex =
ranges.segmentModuleIndices()[pixelModuleIndex] + i_pLS;
871 if (segmentsPixel.isDup()[i_pLS])
873 if (segmentsPixel.partOfPT5()[i_pLS])
876 short layer2_adjustment;
877 if (
modules.layers()[tripletLowerModuleIndex] == 1) {
878 layer2_adjustment = 1;
880 else if (
modules.layers()[tripletLowerModuleIndex] == 2) {
881 layer2_adjustment = 0;
888 for (
unsigned int outerTripletArrayIndex = globalThreadIdx[2]; outerTripletArrayIndex < nOuterTriplets;
889 outerTripletArrayIndex += gridThreadExtent[2]) {
890 unsigned int outerTripletIndex =
891 ranges.tripletModuleIndices()[tripletLowerModuleIndex] + outerTripletArrayIndex;
892 if (
modules.moduleType()[triplets.lowerModuleIndices()[outerTripletIndex][1]] ==
TwoS)
895 if (triplets.partOfPT5()[outerTripletIndex])
898 float pixelRadius, tripletRadius, rPhiChiSquared, rzChiSquared, rPhiChiSquaredInwards, centerX, centerY;
914 rPhiChiSquaredInwards);
918 mds.anchorPhi()[segments
919 .mdIndices()[triplets.segmentIndices()[outerTripletIndex][0]][layer2_adjustment]];
921 mds.anchorEta()[segments
922 .mdIndices()[triplets.segmentIndices()[outerTripletIndex][0]][layer2_adjustment]];
923 float eta_pix = segmentsPixel.eta()[i_pLS];
924 float phi_pix = segmentsPixel.phi()[i_pLS];
925 float pt = segmentsPixel.ptIn()[i_pLS];
926 float score = rPhiChiSquared + rPhiChiSquaredInwards;
927 unsigned int totOccupancyPixelTriplets =
928 alpaka::atomicAdd(acc, &pixelTriplets.totOccupancyPixelTriplets(), 1u, alpaka::hierarchy::Threads{});
931 printf(
"Pixel Triplet excess alert!\n");
934 unsigned int pixelTripletIndex =
935 alpaka::atomicAdd(acc, &pixelTriplets.nPixelTriplets(), 1u, alpaka::hierarchy::Threads{});
947 rPhiChiSquaredInwards,
956 triplets.partOfPT3()[outerTripletIndex] =
true;
965 template <
typename TAcc>
976 betaOut += alpaka::math::copysign(
984 if (betaIn * betaOut > 0.
f and
989 const float betaInUpd =
991 alpaka::math::copysign(
996 const float betaOutUpd =
998 alpaka::math::copysign(
1003 betaAv = 0.5f * (betaInUpd + betaOutUpd);
1006 const float pt_beta_inv =
1009 betaIn += alpaka::math::copysign(
1013 betaOut += alpaka::math::copysign(
1018 betaAv = 0.5f * (betaIn + betaOut);
1026 const float betaInUpd =
1028 alpaka::math::copysign(
1033 const float betaOutUpd =
1035 alpaka::math::copysign(
1042 ? (0.5
f * (betaInUpd + betaOutUpd))
1047 betaIn += alpaka::math::copysign(
1052 betaOut += alpaka::math::copysign(
1058 betaAv = 0.5f * (betaIn + betaOut);
1064 template <
typename TAcc>
1071 uint16_t pixelModuleIndex,
1072 uint16_t outerInnerLowerModuleIndex,
1073 uint16_t outerOuterLowerModuleIndex,
1074 unsigned int innerSegmentIndex,
1075 unsigned int outerSegmentIndex,
1076 unsigned int firstMDIndex,
1077 unsigned int secondMDIndex,
1078 unsigned int thirdMDIndex,
1079 unsigned int fourthMDIndex) {
1080 float dPhi, betaIn, betaOut, pt_beta, zLo, zHi, zLoPointed, zHiPointed, dPhiCut, betaOutCut;
1082 bool isPS_OutLo = (
modules.moduleType()[outerInnerLowerModuleIndex] ==
PS);
1084 float rt_InLo = mds.anchorRt()[firstMDIndex];
1085 float rt_InUp = mds.anchorRt()[secondMDIndex];
1086 float rt_OutLo = mds.anchorRt()[thirdMDIndex];
1087 float rt_OutUp = mds.anchorRt()[fourthMDIndex];
1089 float z_InUp = mds.anchorZ()[secondMDIndex];
1090 float z_OutLo = mds.anchorZ()[thirdMDIndex];
1092 float x_InLo = mds.anchorX()[firstMDIndex];
1093 float x_InUp = mds.anchorX()[secondMDIndex];
1094 float x_OutLo = mds.anchorX()[thirdMDIndex];
1095 float x_OutUp = mds.anchorX()[fourthMDIndex];
1097 float y_InLo = mds.anchorY()[firstMDIndex];
1098 float y_InUp = mds.anchorY()[secondMDIndex];
1099 float y_OutLo = mds.anchorY()[thirdMDIndex];
1100 float y_OutUp = mds.anchorY()[fourthMDIndex];
1102 float rt_InOut = rt_InUp;
1107 unsigned int pixelSegmentArrayIndex = innerSegmentIndex -
ranges.segmentModuleIndices()[pixelModuleIndex];
1108 float ptIn = segmentsPixel.ptIn()[pixelSegmentArrayIndex];
1110 float px = segmentsPixel.px()[pixelSegmentArrayIndex];
1111 float py = segmentsPixel.py()[pixelSegmentArrayIndex];
1112 float pz = segmentsPixel.pz()[pixelSegmentArrayIndex];
1113 float ptErr = segmentsPixel.ptErr()[pixelSegmentArrayIndex];
1114 float etaErr = segmentsPixel.etaErr()[pixelSegmentArrayIndex];
1118 float alpha1GeV_OutLo =
1120 const float rtRatio_OutLoInOut =
1121 rt_OutLo / rt_InOut;
1125 const float zpitch_InLo = 0.05f;
1126 const float zpitch_InOut = 0.05f;
1128 float zGeom = zpitch_InLo + zpitch_OutLo;
1129 zHi = z_InUp + (z_InUp +
kDeltaZLum) * (rtRatio_OutLoInOut - 1.
f) * (z_InUp < 0.f ? 1.f : dzDrtScale) +
1130 (zpitch_InOut + zpitch_OutLo);
1131 zLo = z_InUp + (z_InUp -
kDeltaZLum) * (rtRatio_OutLoInOut - 1.
f) * (z_InUp > 0.f ? 1.f : dzDrtScale) -
1132 (zpitch_InOut + zpitch_OutLo);
1134 if ((z_OutLo < zLo) || (z_OutLo > zHi))
1137 const float cosh2Eta = 1.f + (pz * pz) / (ptIn * ptIn);
1139 const float drt_OutLo_InUp = (rt_OutLo - rt_InUp);
1143 float drt_InSeg = rt_InOut - rt_InLo;
1145 const float thetaMuls2 =
1147 const float muls2 = thetaMuls2 * 9.f / (
ptCut *
ptCut) * 16.
f;
1149 float dzErr = (drt_OutLo_InUp * drt_OutLo_InUp) * (etaErr * etaErr) * cosh2Eta;
1150 dzErr += 0.03f * 0.03f;
1152 dzErr += muls2 * (drt_OutLo_InUp * drt_OutLo_InUp) / 3.
f * cosh2Eta;
1153 dzErr += zGeom * zGeom;
1156 const float dzDrIn = pz / ptIn;
1157 const float zWindow =
dzErr / drt_InSeg * drt_OutLo_InUp + zGeom;
1158 const float dzMean = dzDrIn * drt_OutLo_InUp *
1162 zLoPointed = z_InUp + dzMean - zWindow;
1163 zHiPointed = z_InUp + dzMean + zWindow;
1165 if ((z_OutLo < zLoPointed) || (z_OutLo > zHiPointed))
1168 const float pvOffset = 0.1f / rt_OutLo;
1172 float midPointX = 0.5f * (x_InLo + x_OutLo);
1173 float midPointY = 0.5f * (y_InLo + y_OutLo);
1175 float diffX = x_OutLo - x_InLo;
1176 float diffY = y_OutLo - y_InLo;
1178 dPhi =
deltaPhi(acc, midPointX, midPointY, diffX, diffY);
1185 float alpha_InLo =
__H2F(segments.dPhiChanges()[innerSegmentIndex]);
1186 float alpha_OutLo =
__H2F(segments.dPhiChanges()[outerSegmentIndex]);
1188 bool isEC_lastLayer =
modules.subdets()[outerOuterLowerModuleIndex] ==
Endcap and
1189 modules.moduleType()[outerOuterLowerModuleIndex] ==
TwoS;
1191 float alpha_OutUp, alpha_OutUp_highEdge, alpha_OutUp_lowEdge;
1192 alpha_OutUp =
deltaPhi(acc, x_OutUp, y_OutUp, x_OutUp - x_OutLo, y_OutUp - y_OutLo);
1194 alpha_OutUp_highEdge = alpha_OutUp;
1195 alpha_OutUp_lowEdge = alpha_OutUp;
1197 float tl_axis_x = x_OutUp - x_InUp;
1198 float tl_axis_y = y_OutUp - y_InUp;
1200 float tl_axis_highEdge_x = tl_axis_x;
1201 float tl_axis_highEdge_y = tl_axis_y;
1203 float tl_axis_lowEdge_x = tl_axis_x;
1204 float tl_axis_lowEdge_y = tl_axis_y;
1207 float betaInRHmin = betaIn;
1208 float betaInRHmax = betaIn;
1210 betaOut = -alpha_OutUp +
deltaPhi(acc, x_OutUp, y_OutUp, tl_axis_x, tl_axis_y);
1212 float betaOutRHmin = betaOut;
1213 float betaOutRHmax = betaOut;
1215 if (isEC_lastLayer) {
1216 alpha_OutUp_highEdge =
deltaPhi(acc,
1217 mds.anchorHighEdgeX()[fourthMDIndex],
1218 mds.anchorHighEdgeY()[fourthMDIndex],
1219 mds.anchorHighEdgeX()[fourthMDIndex] - x_OutLo,
1220 mds.anchorHighEdgeY()[fourthMDIndex] - y_OutLo);
1221 alpha_OutUp_lowEdge =
deltaPhi(acc,
1222 mds.anchorLowEdgeX()[fourthMDIndex],
1223 mds.anchorLowEdgeY()[fourthMDIndex],
1224 mds.anchorLowEdgeX()[fourthMDIndex] - x_OutLo,
1225 mds.anchorLowEdgeY()[fourthMDIndex] - y_OutLo);
1227 tl_axis_highEdge_x = mds.anchorHighEdgeX()[fourthMDIndex] - x_InUp;
1228 tl_axis_highEdge_y = mds.anchorHighEdgeY()[fourthMDIndex] - y_InUp;
1229 tl_axis_lowEdge_x = mds.anchorLowEdgeX()[fourthMDIndex] - x_InUp;
1230 tl_axis_lowEdge_y = mds.anchorLowEdgeY()[fourthMDIndex] - y_InUp;
1232 betaOutRHmin = -alpha_OutUp_highEdge +
deltaPhi(acc,
1233 mds.anchorHighEdgeX()[fourthMDIndex],
1234 mds.anchorHighEdgeY()[fourthMDIndex],
1236 tl_axis_highEdge_y);
1237 betaOutRHmax = -alpha_OutUp_lowEdge +
deltaPhi(acc,
1238 mds.anchorLowEdgeX()[fourthMDIndex],
1239 mds.anchorLowEdgeY()[fourthMDIndex],
1245 float drt_tl_axis =
alpaka::math::sqrt(acc, tl_axis_x * tl_axis_x + tl_axis_y * tl_axis_y);
1248 const float rt_InSeg =
1249 alpaka::math::sqrt(acc, (x_InUp - x_InLo) * (x_InUp - x_InLo) + (y_InUp - y_InLo) * (y_InUp - y_InLo));
1252 float betaAv = 0.5f * (betaIn + betaOut);
1256 int lOut = isEC_lastLayer ? 11 : 5;
1258 alpaka::math::sqrt(acc, (x_OutUp - x_OutLo) * (x_OutUp - x_OutLo) + (y_OutUp - y_OutLo) * (y_OutUp - y_OutLo));
1259 float sdOut_d = rt_OutUp - rt_OutLo;
1263 const float betaInMMSF = (
alpaka::math::abs(acc, betaInRHmin + betaInRHmax) > 0)
1266 const float betaOutMMSF = (
alpaka::math::abs(acc, betaOutRHmin + betaOutRHmax) > 0)
1269 betaInRHmin *= betaInMMSF;
1270 betaInRHmax *= betaInMMSF;
1271 betaOutRHmin *= betaOutMMSF;
1272 betaOutRHmax *= betaOutMMSF;
1276 const float dBetaMuls2 = thetaMuls2 * 16.f / (min_ptBeta_ptBetaMax * min_ptBeta_ptBetaMax);
1277 const float alphaInAbsReg =
1281 const float alphaOutAbsReg =
1287 const float dBetaLum2 = (dBetaInLum + dBetaOutLum) * (dBetaInLum + dBetaOutLum);
1290 const float dBetaRIn2 = 0;
1292 float dBetaROut = 0;
1293 if (isEC_lastLayer) {
1295 mds.anchorHighEdgeX()[fourthMDIndex] * mds.anchorHighEdgeX()[fourthMDIndex] +
1296 mds.anchorHighEdgeY()[fourthMDIndex] * mds.anchorHighEdgeY()[fourthMDIndex]) -
1298 mds.anchorLowEdgeX()[fourthMDIndex] * mds.anchorLowEdgeX()[fourthMDIndex] +
1299 mds.anchorLowEdgeY()[fourthMDIndex] * mds.anchorLowEdgeY()[fourthMDIndex])) *
1300 sinDPhi / drt_tl_axis;
1303 const float dBetaROut2 = dBetaROut * dBetaROut;
1313 const float dBetaCut2 =
1314 (dBetaRes * dBetaRes * 2.0f + dBetaMuls2 + dBetaLum2 + dBetaRIn2 + dBetaROut2 +
1318 float dBeta = betaIn - betaOut;
1319 return dBeta * dBeta <= dBetaCut2;
1322 template <
typename TAcc>
1329 uint16_t pixelModuleIndex,
1330 uint16_t outerInnerLowerModuleIndex,
1331 uint16_t outerOuterLowerModuleIndex,
1332 unsigned int innerSegmentIndex,
1333 unsigned int outerSegmentIndex,
1334 unsigned int firstMDIndex,
1335 unsigned int secondMDIndex,
1336 unsigned int thirdMDIndex,
1337 unsigned int fourthMDIndex) {
1338 float dPhi, betaIn, betaOut, pt_beta, rtLo, rtHi, dPhiCut, betaOutCut;
1340 bool isPS_OutLo = (
modules.moduleType()[outerInnerLowerModuleIndex] ==
PS);
1342 float z_InUp = mds.anchorZ()[secondMDIndex];
1343 float z_OutLo = mds.anchorZ()[thirdMDIndex];
1345 if (z_InUp * z_OutLo <= 0)
1348 float rt_InLo = mds.anchorRt()[firstMDIndex];
1349 float rt_InUp = mds.anchorRt()[secondMDIndex];
1350 float rt_OutLo = mds.anchorRt()[thirdMDIndex];
1351 float rt_OutUp = mds.anchorRt()[fourthMDIndex];
1353 float x_InLo = mds.anchorX()[firstMDIndex];
1354 float x_InUp = mds.anchorX()[secondMDIndex];
1355 float x_OutLo = mds.anchorX()[thirdMDIndex];
1356 float x_OutUp = mds.anchorX()[fourthMDIndex];
1358 float y_InLo = mds.anchorY()[firstMDIndex];
1359 float y_InUp = mds.anchorY()[secondMDIndex];
1360 float y_OutLo = mds.anchorY()[thirdMDIndex];
1361 float y_OutUp = mds.anchorY()[fourthMDIndex];
1363 unsigned int pixelSegmentArrayIndex = innerSegmentIndex -
ranges.segmentModuleIndices()[pixelModuleIndex];
1365 float ptIn = segmentsPixel.ptIn()[pixelSegmentArrayIndex];
1367 float px = segmentsPixel.px()[pixelSegmentArrayIndex];
1368 float py = segmentsPixel.py()[pixelSegmentArrayIndex];
1369 float pz = segmentsPixel.pz()[pixelSegmentArrayIndex];
1370 float ptErr = segmentsPixel.ptErr()[pixelSegmentArrayIndex];
1371 float etaErr = segmentsPixel.etaErr()[pixelSegmentArrayIndex];
1376 const float zpitch_InLo = 0.05f;
1378 float zGeom = zpitch_InLo + zpitch_OutLo;
1383 const float dLum = alpaka::math::copysign(acc,
kDeltaZLum, z_InUp);
1384 bool isOutSgInnerMDPS =
modules.moduleType()[outerInnerLowerModuleIndex] ==
PS;
1386 const float rtGeom1 = isOutSgInnerMDPS
1389 const float zGeom1 = alpaka::math::copysign(acc, zGeom, z_InUp);
1390 rtLo = rt_InUp * (1.f + (z_OutLo - z_InUp - zGeom1) / (z_InUp + zGeom1 + dLum) / dzDrtScale) -
1393 float zInForHi = z_InUp - zGeom1 - dLum;
1394 if (zInForHi * z_InUp < 0)
1395 zInForHi = alpaka::math::copysign(acc, 0.1
f, z_InUp);
1396 rtHi = rt_InUp * (1.f + (z_OutLo - z_InUp + zGeom1) / zInForHi) + rtGeom1;
1399 if ((rt_OutLo < rtLo) || (rt_OutLo > rtHi))
1403 const float cosh2Eta = 1.f + (pz * pz) / (ptIn * ptIn);
1404 const float multDzDr2 = (dzOutInAbs * dzOutInAbs) * cosh2Eta / ((cosh2Eta - 1.
f) * (cosh2Eta - 1.f));
1406 const float thetaMuls2 =
1408 const float muls2 = thetaMuls2 * 9.f / (
ptCut *
ptCut) * 16.
f;
1410 float drtErr = (etaErr * etaErr) * multDzDr2;
1411 drtErr += 0.03f * 0.03f;
1413 drtErr += muls2 * multDzDr2 / 3.f * cosh2Eta;
1417 const float drt_OutLo_InUp = (rt_OutLo - rt_InUp);
1419 const float rtWindow = drtErr + rtGeom1;
1420 const float drtMean = drtDzIn * dzOutInAbs *
1423 const float rtLo_point = rt_InUp + drtMean - rtWindow;
1424 const float rtHi_point = rt_InUp + drtMean + rtWindow;
1427 if ((rt_OutLo < rtLo_point) || (rt_OutLo > rtHi_point))
1430 const float alpha1GeV_OutLo =
1432 const float pvOffset = 0.1f / rt_OutLo;
1435 float midPointX = 0.5f * (x_InLo + x_OutLo);
1436 float midPointY = 0.5f * (y_InLo + y_OutLo);
1438 float diffX = x_OutLo - x_InLo;
1439 float diffY = y_OutLo - y_InLo;
1441 dPhi =
deltaPhi(acc, midPointX, midPointY, diffX, diffY);
1447 float alpha_InLo =
__H2F(segments.dPhiChanges()[innerSegmentIndex]);
1448 float alpha_OutLo =
__H2F(segments.dPhiChanges()[outerSegmentIndex]);
1450 bool isEC_lastLayer =
modules.subdets()[outerOuterLowerModuleIndex] ==
Endcap and
1451 modules.moduleType()[outerOuterLowerModuleIndex] ==
TwoS;
1453 float alpha_OutUp, alpha_OutUp_highEdge, alpha_OutUp_lowEdge;
1455 alpha_OutUp =
deltaPhi(acc, x_OutUp, y_OutUp, x_OutUp - x_OutLo, y_OutUp - y_OutLo);
1456 alpha_OutUp_highEdge = alpha_OutUp;
1457 alpha_OutUp_lowEdge = alpha_OutUp;
1459 float tl_axis_x = x_OutUp - x_InUp;
1460 float tl_axis_y = y_OutUp - y_InUp;
1462 float tl_axis_highEdge_x = tl_axis_x;
1463 float tl_axis_highEdge_y = tl_axis_y;
1465 float tl_axis_lowEdge_x = tl_axis_x;
1466 float tl_axis_lowEdge_y = tl_axis_y;
1469 float betaInRHmin = betaIn;
1470 float betaInRHmax = betaIn;
1472 betaOut = -alpha_OutUp +
deltaPhi(acc, x_OutUp, y_OutUp, tl_axis_x, tl_axis_y);
1473 float betaOutRHmin = betaOut;
1474 float betaOutRHmax = betaOut;
1476 if (isEC_lastLayer) {
1477 alpha_OutUp_highEdge =
deltaPhi(acc,
1478 mds.anchorHighEdgeX()[fourthMDIndex],
1479 mds.anchorHighEdgeY()[fourthMDIndex],
1480 mds.anchorHighEdgeX()[fourthMDIndex] - x_OutLo,
1481 mds.anchorHighEdgeY()[fourthMDIndex] - y_OutLo);
1482 alpha_OutUp_lowEdge =
deltaPhi(acc,
1483 mds.anchorLowEdgeX()[fourthMDIndex],
1484 mds.anchorLowEdgeY()[fourthMDIndex],
1485 mds.anchorLowEdgeX()[fourthMDIndex] - x_OutLo,
1486 mds.anchorLowEdgeY()[fourthMDIndex] - y_OutLo);
1488 tl_axis_highEdge_x = mds.anchorHighEdgeX()[fourthMDIndex] - x_InUp;
1489 tl_axis_highEdge_y = mds.anchorHighEdgeY()[fourthMDIndex] - y_InUp;
1490 tl_axis_lowEdge_x = mds.anchorLowEdgeX()[fourthMDIndex] - x_InUp;
1491 tl_axis_lowEdge_y = mds.anchorLowEdgeY()[fourthMDIndex] - y_InUp;
1493 betaOutRHmin = -alpha_OutUp_highEdge +
deltaPhi(acc,
1494 mds.anchorHighEdgeX()[fourthMDIndex],
1495 mds.anchorHighEdgeY()[fourthMDIndex],
1497 tl_axis_highEdge_y);
1498 betaOutRHmax = -alpha_OutUp_lowEdge +
deltaPhi(acc,
1499 mds.anchorLowEdgeX()[fourthMDIndex],
1500 mds.anchorLowEdgeY()[fourthMDIndex],
1506 float drt_tl_axis =
alpaka::math::sqrt(acc, tl_axis_x * tl_axis_x + tl_axis_y * tl_axis_y);
1508 const float rt_InSeg =
1509 alpaka::math::sqrt(acc, (x_InUp - x_InLo) * (x_InUp - x_InLo) + (y_InUp - y_InLo) * (y_InUp - y_InLo));
1511 float betaAv = 0.5f * (betaIn + betaOut);
1515 int lOut = isEC_lastLayer ? 11 : 5;
1517 alpaka::math::sqrt(acc, (x_OutUp - x_OutLo) * (x_OutUp - x_OutLo) + (y_OutUp - y_OutLo) * (y_OutUp - y_OutLo));
1518 float sdOut_d = rt_OutUp - rt_OutLo;
1522 const float betaInMMSF = (
alpaka::math::abs(acc, betaInRHmin + betaInRHmax) > 0)
1525 const float betaOutMMSF = (
alpaka::math::abs(acc, betaOutRHmin + betaOutRHmax) > 0)
1528 betaInRHmin *= betaInMMSF;
1529 betaInRHmax *= betaInMMSF;
1530 betaOutRHmin *= betaOutMMSF;
1531 betaOutRHmax *= betaOutMMSF;
1535 const float dBetaMuls2 = thetaMuls2 * 16.f / (min_ptBeta_ptBetaMax * min_ptBeta_ptBetaMax);
1537 const float alphaInAbsReg =
1541 const float alphaOutAbsReg =
1547 const float dBetaLum2 = (dBetaInLum + dBetaOutLum) * (dBetaInLum + dBetaOutLum);
1550 const float dBetaRIn2 = 0;
1552 float dBetaROut = 0;
1553 if (isEC_lastLayer) {
1555 mds.anchorHighEdgeX()[fourthMDIndex] * mds.anchorHighEdgeX()[fourthMDIndex] +
1556 mds.anchorHighEdgeY()[fourthMDIndex] * mds.anchorHighEdgeY()[fourthMDIndex]) -
1558 mds.anchorLowEdgeX()[fourthMDIndex] * mds.anchorLowEdgeX()[fourthMDIndex] +
1559 mds.anchorLowEdgeY()[fourthMDIndex] * mds.anchorLowEdgeY()[fourthMDIndex])) *
1560 sinDPhi / drt_tl_axis;
1563 const float dBetaROut2 = dBetaROut * dBetaROut;
1574 float drt_InSeg = rt_InUp - rt_InLo;
1577 const float dBetaCut2 =
1578 (dBetaRes * dBetaRes * 2.0f + dBetaMuls2 + dBetaLum2 + dBetaRIn2 + dBetaROut2 +
1582 float dBeta = betaIn - betaOut;
1583 return dBeta * dBeta <= dBetaCut2;
ALPAKA_FN_ACC ALPAKA_FN_INLINE bool passRadiusCriterionBBE(TAcc const &acc, float pixelRadius, float pixelRadiusError, float tripletRadius)
ALPAKA_FN_ACC ALPAKA_FN_INLINE bool passRadiusCriterionEEE(TAcc const &acc, float pixelRadius, float pixelRadiusError, float tripletRadius)
ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kR1GeVf
ALPAKA_FN_ACC ALPAKA_FN_INLINE bool passRadiusCriterionBEE(TAcc const &acc, float pixelRadius, float pixelRadiusError, float tripletRadius)
ALPAKA_FN_ACC ALPAKA_FN_INLINE float computeChiSquaredpT3(TAcc const &acc, unsigned int nPoints, float *xs, float *ys, float *delta1, float *delta2, float *slopes, bool *isFlat, float g, float f, float radius)
ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float k2Rinv1GeVf
ALPAKA_FN_ACC ALPAKA_FN_INLINE float computePT3RZChiSquared(TAcc const &acc, ModulesConst modules, const uint16_t *lowerModuleIndices, const float *rtPix, const float *xPix, const float *yPix, const float *zPix, const float *rts, const float *xs, const float *ys, const float *zs, float pixelSegmentPt, float pixelSegmentPx, float pixelSegmentPy, float pixelSegmentPz, int pixelSegmentCharge)
ALPAKA_FN_ACC ALPAKA_FN_INLINE float computePT3RPhiChiSquared(TAcc const &acc, ModulesConst modules, uint16_t *lowerModuleIndices, float g, float f, float radius, float *xs, float *ys)
static const double slope[3]
Sin< T >::type sin(const T &t)
ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kWidthPS
ALPAKA_FN_ACC ALPAKA_FN_INLINE bool runTripletDefaultAlgoPPEE(TAcc const &acc, ModulesConst modules, ObjectRangesConst ranges, MiniDoubletsConst mds, SegmentsConst segments, SegmentsPixelConst segmentsPixel, uint16_t pixelModuleIndex, uint16_t outerInnerLowerModuleIndex, uint16_t outerOuterLowerModuleIndex, unsigned int innerSegmentIndex, unsigned int outerSegmentIndex, unsigned int firstMDIndex, unsigned int secondMDIndex, unsigned int thirdMDIndex, unsigned int fourthMDIndex)
ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kWidth2S
ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kVerticalModuleSlope
constexpr unsigned int n_max_pixel_triplets
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
ALPAKA_FN_ACC ALPAKA_FN_INLINE float computePT3RPhiChiSquaredInwards(float g, float f, float r, float *xPix, float *yPix)
ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float ptCut
ALPAKA_FN_ACC ALPAKA_FN_INLINE bool runTripletDefaultAlgoPPBB(TAcc const &acc, ModulesConst modules, ObjectRangesConst ranges, MiniDoubletsConst mds, SegmentsConst segments, SegmentsPixelConst segmentsPixel, uint16_t pixelModuleIndex, uint16_t outerInnerLowerModuleIndex, uint16_t outerOuterLowerModuleIndex, unsigned int innerSegmentIndex, unsigned int outerSegmentIndex, unsigned int firstMDIndex, unsigned int secondMDIndex, unsigned int thirdMDIndex, unsigned int fourthMDIndex)
ALPAKA_FN_ACC ALPAKA_FN_INLINE bool runPixelTripletDefaultAlgo(TAcc const &acc, ModulesConst modules, ObjectRangesConst ranges, MiniDoubletsConst mds, SegmentsConst segments, SegmentsPixelConst segmentsPixel, TripletsConst triplets, unsigned int pixelSegmentIndex, unsigned int tripletIndex, float &pixelRadius, float &tripletRadius, float ¢erX, float ¢erY, float &rzChiSquared, float &rPhiChiSquared, float &rPhiChiSquaredInwards, bool runChiSquaredCuts=true)
TripletsSoA::ConstView TripletsConst
ALPAKA_FN_ACC void operator()(TAcc const &acc, ModulesConst modules, ModulesPixelConst modulesPixel, ObjectRangesConst ranges, MiniDoubletsConst mds, SegmentsConst segments, SegmentsPixelConst segmentsPixel, Triplets triplets, TripletsOccupancyConst tripletsOccupancy, PixelTriplets pixelTriplets, unsigned int *connectedPixelSize, unsigned int *connectedPixelIndex, unsigned int nPixelSegments) const
ALPAKA_FN_ACC ALPAKA_FN_INLINE bool passPT3RPhiChiSquaredCuts(ModulesConst modules, uint16_t lowerModuleIndex1, uint16_t lowerModuleIndex2, uint16_t lowerModuleIndex3, float chiSquared)
ALPAKA_FN_HOST_ACC ALPAKA_FN_INLINE float phi(TAcc const &acc, float x, float y)
ModulesPixelSoA::ConstView ModulesPixelConst
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
Cos< T >::type cos(const T &t)
Tan< T >::type tan(const T &t)
Abs< T >::type abs(const T &t)
ALPAKA_FN_ACC ALPAKA_FN_INLINE bool passPT3RZChiSquaredCuts(ModulesConst modules, uint16_t lowerModuleIndex1, uint16_t lowerModuleIndex2, uint16_t lowerModuleIndex3, float rzChiSquared)
MiniDoubletsSoA::ConstView MiniDoubletsConst
ModulesSoA::ConstView ModulesConst
PixelTripletsSoA::View PixelTriplets
SegmentsPixelSoA::ConstView SegmentsPixelConst
ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kMulsInGeV
ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kPixelPSZpitch
ALPAKA_FN_ACC ALPAKA_FN_INLINE bool checkIntervalOverlappT3(float firstMin, float firstMax, float secondMin, float secondMax)
ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kPt_betaMax
ALPAKA_FN_ACC ALPAKA_FN_INLINE bool passRadiusCriterion(TAcc const &acc, ModulesConst modules, float pixelRadius, float pixelRadiusError, float tripletRadius, int16_t lowerModuleIndex, uint16_t middleModuleIndex, uint16_t upperModuleIndex)
ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kPi
ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kStrip2SZpitch
TripletsSoA::View Triplets
ALPAKA_FN_ACC ALPAKA_FN_INLINE bool passPT3RPhiChiSquaredInwardsCuts(ModulesConst modules, uint16_t lowerModuleIndex1, uint16_t lowerModuleIndex2, uint16_t lowerModuleIndex3, float chiSquared)
ALPAKA_FN_ACC ALPAKA_FN_INLINE bool passRadiusCriterionBBB(TAcc const &acc, float pixelRadius, float pixelRadiusError, float tripletRadius)
TripletsOccupancySoA::ConstView TripletsOccupancyConst
ObjectRangesSoA::ConstView ObjectRangesConst
ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kSinAlphaMax
SegmentsSoA::ConstView SegmentsConst
ALPAKA_STATIC_ACC_MEM_GLOBAL constexpr float kDeltaZLum
ALPAKA_FN_ACC ALPAKA_FN_INLINE bool runPixelTrackletDefaultAlgopT3(TAcc const &acc, ModulesConst modules, ObjectRangesConst ranges, MiniDoubletsConst mds, SegmentsConst segments, SegmentsPixelConst segmentsPixel, uint16_t pixelLowerModuleIndex, uint16_t outerInnerLowerModuleIndex, uint16_t outerOuterLowerModuleIndex, unsigned int innerSegmentIndex, unsigned int outerSegmentIndex)
ALPAKA_FN_HOST_ACC ALPAKA_FN_INLINE float deltaPhi(TAcc const &acc, float x1, float y1, float x2, float y2)
T1 atomicAdd(T1 *a, T2 b)
ALPAKA_FN_ACC ALPAKA_FN_INLINE void addPixelTripletToMemory(MiniDoubletsConst mds, SegmentsConst segments, TripletsConst triplets, PixelTriplets pixelTriplets, unsigned int pixelSegmentIndex, unsigned int tripletIndex, float pixelRadius, float tripletRadius, float centerX, float centerY, float rPhiChiSquared, float rPhiChiSquaredInwards, float rzChiSquared, unsigned int pixelTripletIndex, float pt, float eta, float phi, float eta_pix, float phi_pix, float score)
ALPAKA_FN_ACC ALPAKA_FN_INLINE void runDeltaBetaIterationspT3(TAcc const &acc, float &betaIn, float &betaOut, float betaAv, float &pt_beta, float sdIn_dr, float sdOut_dr, float dr, float lIn)
ALPAKA_FN_HOST_ACC ALPAKA_FN_INLINE float eta(TAcc const &acc, float x, float y, float z)