Merge pull request #40 from YonsiG/CMSSW_14_1_0_pre3_LST_X_LSTCore_realfiles_TrackLooperPR410

LST object build upgrade: change pT5 rz chi2 cut definition into helix approximation

Author: VourMa

RecoTracker/LSTCore/src/alpaka/PixelTriplet.h

@@ -755,8 +755,7 @@ namespace SDL {
     float z1 = zPix[1] / 100;
     float r1 = rtPix[1] / 100;
 
-    float Bz = SDL::magnetic_field;
-    float a = -0.299792 * Bz * charge;
+    float a = -100 / SDL::kR1GeVf * charge;
 
     for (size_t i = 0; i < 3; i++) {
       float zsi = zs[i] / 100;
@@ -776,6 +775,7 @@ namespace SDL {
         float x = x1 + Px / a * alpaka::math::sin(acc, rou * s) - Py / a * (1 - alpaka::math::cos(acc, rou * s));
         float y = y1 + Py / a * alpaka::math::sin(acc, rou * s) + Px / a * (1 - alpaka::math::cos(acc, rou * s));
         diffr = alpaka::math::abs(acc, rtsi - alpaka::math::sqrt(acc, x * x + y * y)) * 100;
+        residual = diffr;
       }
 
       if (moduleSubdet == SDL::Barrel) {
@@ -792,10 +792,9 @@ namespace SDL {
         float diffz1 = alpaka::math::abs(acc, solz1 - zsi) * 100;
         float diffz2 = alpaka::math::abs(acc, solz2 - zsi) * 100;
         diffz = alpaka::math::min(acc, diffz1, diffz2);
+        residual = diffz;
       }
 
-      residual = moduleSubdet == SDL::Barrel ? diffz : diffr;
-
       //PS Modules
       if (moduleType == 0) {
         error = 0.15f;
@@ -2053,63 +2052,61 @@ namespace SDL {
                        5 * (modulesInGPU.subdets[lowerModuleIndex5] == SDL::Endcap and
                             modulesInGPU.moduleType[lowerModuleIndex5] == SDL::TwoS);
 
+    // This slides shows the cut threshold definition. The comments below in the code, e.g, "cat 10", is consistent with the region separation in the slides
+    // https://indico.cern.ch/event/1410985/contributions/5931017/attachments/2875400/5035406/helix%20approxi%20for%20pT5%20rzchi2%20new%20results%20versions.pdf
     if (layer1 == 1 and layer2 == 2 and layer3 == 3) {
-      if (layer4 == 12 and layer5 == 13) {
-        return rzChiSquared < 451.141f;
-      } else if (layer4 == 4 and layer5 == 12) {
-        return rzChiSquared < 392.654f;
-      } else if (layer4 == 4 and layer5 == 5) {
-        return rzChiSquared < 225.322f;
-      } else if (layer4 == 7 and layer5 == 13) {
-        return rzChiSquared < 595.546f;
-      } else if (layer4 == 7 and layer5 == 8) {
-        return rzChiSquared < 196.111f;
+      if (layer4 == 12 and layer5 == 13) {  // cat 10
+        return rzChiSquared < 14.031f;
+      } else if (layer4 == 4 and layer5 == 12) {  // cat 12
+        return rzChiSquared < 8.760f;
+      } else if (layer4 == 4 and layer5 == 5) {  // cat 11
+        return rzChiSquared < 3.607f;
+      } else if (layer4 == 7 and layer5 == 13) {  // cat 9
+        return rzChiSquared < 16.620;
+      } else if (layer4 == 7 and layer5 == 8) {  // cat 8
+        return rzChiSquared < 17.910f;
       }
     } else if (layer1 == 1 and layer2 == 2 and layer3 == 7) {
-      if (layer4 == 13 and layer5 == 14) {
-        return rzChiSquared < 297.446f;
-      } else if (layer4 == 8 and layer5 == 14) {
-        return rzChiSquared < 451.141f;
-      } else if (layer4 == 8 and layer5 == 9) {
-        return rzChiSquared < 518.339f;
+      if (layer4 == 13 and layer5 == 14) {  // cat 7
+        return rzChiSquared < 8.950f;
+      } else if (layer4 == 8 and layer5 == 14) {  // cat 6
+        return rzChiSquared < 14.837f;
+      } else if (layer4 == 8 and layer5 == 9) {  // cat 5
+        return rzChiSquared < 18.519f;
       }
     } else if (layer1 == 1 and layer2 == 7 and layer3 == 8) {
-      if (layer4 == 9 and layer5 == 10) {
-        return rzChiSquared < 341.75f;
-      } else if (layer4 == 9 and layer5 == 15) {
-        return rzChiSquared < 341.75f;
+      if (layer4 == 9 and layer5 == 10) {  // cat 3
+        return rzChiSquared < 15.093f;
+      } else if (layer4 == 9 and layer5 == 15) {  // cat 4
+        return rzChiSquared < 11.200f;
       }
     } else if (layer1 == 2 and layer2 == 3 and layer3 == 4) {
-      if (layer4 == 12 and layer5 == 13) {
-        return rzChiSquared < 392.655f;
-      } else if (layer4 == 5 and layer5 == 12) {
-        return rzChiSquared < 341.75f;
-      } else if (layer4 == 5 and layer5 == 6) {
-        return rzChiSquared < 112.537f;
+      if (layer4 == 12 and layer5 == 13) {  // cat 20
+        return rzChiSquared < 12.868f;
+      } else if (layer4 == 5 and layer5 == 12) {  // cat 19
+        return rzChiSquared < 6.128f;
+      } else if (layer4 == 5 and layer5 == 6) {  // cat 18
+        return rzChiSquared < 2.987f;
       }
     } else if (layer1 == 2 and layer2 == 3 and layer4 == 7) {
-      if (layer4 == 13 and layer5 == 14) {
-        return rzChiSquared < 595.545f;
-      } else if (layer4 == 8 and layer5 == 14) {
-        return rzChiSquared < 74.198f;
+      if (layer4 == 13 and layer5 == 14) {  // cat 17
+        return rzChiSquared < 19.446f;
+      } else if (layer4 == 8 and layer5 == 14) {  // cat 16
+        return rzChiSquared < 17.520f;
       }
     } else if (layer1 == 2 and layer2 == 7 and layer3 == 8) {
-      if (layer4 == 14 and layer5 == 15) {
-        return rzChiSquared < 518.339f;
-      } else if (layer4 == 9 and layer5 == 10) {
-        return rzChiSquared < 8.046f;
-      } else if (layer4 == 9 and layer5 == 15) {
-        return rzChiSquared < 451.141f;
+      if (layer4 == 14 and layer5 == 15) {  // cat 15
+        return rzChiSquared < 14.71f;
+      } else if (layer4 == 9 and layer5 == 15) {  // cat 14
+        return rzChiSquared < 18.213f;
       }
-    } else if (layer1 == 3 and layer2 == 7 and layer3 == 8 and layer4 == 14 and layer5 == 15) {
-      return rzChiSquared < 56.207f;
     } else if (layer1 == 7 and layer2 == 8 and layer3 == 9) {
-      if (layer4 == 10 and layer5 == 11) {
-        return rzChiSquared < 64.578f;
-      } else if (layer4 == 10 and layer5 == 16) {
-        return rzChiSquared < 85.250f;
-      } else if (layer4 == 15 and layer5 == 16) {
-        return rzChiSquared < 85.250f;
+      if (layer4 == 10 and layer5 == 11) {  // cat 0
+        return rzChiSquared < 10.016f;
+      } else if (layer4 == 10 and layer5 == 16) {  // cat 1
+        return rzChiSquared < 87.671f;
+      } else if (layer4 == 15 and layer5 == 16) {  // cat 2
+        return rzChiSquared < 5.844f;
       }
     }
     return true;
@@ -2534,8 +2531,11 @@ namespace SDL {
     uint16_t lowerModuleIndices[5] = {
         lowerModuleIndex1, lowerModuleIndex2, lowerModuleIndex3, lowerModuleIndex4, lowerModuleIndex5};
 
-    float zPix[2] = {mdsInGPU.anchorZ[pixelInnerMDIndex], mdsInGPU.anchorZ[pixelOuterMDIndex]};
     float rtPix[2] = {mdsInGPU.anchorRt[pixelInnerMDIndex], mdsInGPU.anchorRt[pixelOuterMDIndex]};
+    float xPix[2] = {mdsInGPU.anchorX[pixelInnerMDIndex], mdsInGPU.anchorX[pixelOuterMDIndex]};
+    float yPix[2] = {mdsInGPU.anchorY[pixelInnerMDIndex], mdsInGPU.anchorY[pixelOuterMDIndex]};
+    float zPix[2] = {mdsInGPU.anchorZ[pixelInnerMDIndex], mdsInGPU.anchorZ[pixelOuterMDIndex]};
+
     float zs[5] = {mdsInGPU.anchorZ[firstMDIndex],
                    mdsInGPU.anchorZ[secondMDIndex],
                    mdsInGPU.anchorZ[thirdMDIndex],
@@ -2547,9 +2547,32 @@ namespace SDL {
                     mdsInGPU.anchorRt[fourthMDIndex],
                     mdsInGPU.anchorRt[fifthMDIndex]};
 
-    rzChiSquared = computePT5RZChiSquared(acc, modulesInGPU, lowerModuleIndices, rtPix, zPix, rts, zs);
+    float pixelSegmentPt = segmentsInGPU.ptIn[pixelSegmentArrayIndex];
+    float pixelSegmentPx = segmentsInGPU.px[pixelSegmentArrayIndex];
+    float pixelSegmentPy = segmentsInGPU.py[pixelSegmentArrayIndex];
+    float pixelSegmentPz = segmentsInGPU.pz[pixelSegmentArrayIndex];
+    int pixelSegmentCharge = segmentsInGPU.charge[pixelSegmentArrayIndex];
+
+    rzChiSquared = 0;
 
-    if (/*pixelRadius*/ 0 < 5.0f * kR1GeVf) {  // FIXME: pixelRadius is not defined yet
+    //get the appropriate centers
+    pixelRadius = segmentsInGPU.circleRadius[pixelSegmentArrayIndex];
+
+    if (pixelRadius < 5.0f * kR1GeVf) {  //only apply r-z chi2 cuts for <5GeV tracks
+      rzChiSquared = computePT5RZChiSquared(acc,
+                                            modulesInGPU,
+                                            lowerModuleIndices,
+                                            rtPix,
+                                            xPix,
+                                            yPix,
+                                            zPix,
+                                            rts,
+                                            zs,
+                                            pixelSegmentPt,
+                                            pixelSegmentPx,
+                                            pixelSegmentPy,
+                                            pixelSegmentPz,
+                                            pixelSegmentCharge);
       pass = pass and passPT5RZChiSquaredCuts(modulesInGPU,
                                               lowerModuleIndex1,
                                               lowerModuleIndex2,
@@ -2573,10 +2596,9 @@ namespace SDL {
                    mdsInGPU.anchorY[fourthMDIndex],
                    mdsInGPU.anchorY[fifthMDIndex]};
 
-    //get the appropriate radii and centers
+    //get the appropriate centers
     centerX = segmentsInGPU.circleCenterX[pixelSegmentArrayIndex];
     centerY = segmentsInGPU.circleCenterY[pixelSegmentArrayIndex];
-    pixelRadius = segmentsInGPU.circleRadius[pixelSegmentArrayIndex];
 
     float T5CenterX = quintupletsInGPU.regressionG[quintupletIndex];
     float T5CenterY = quintupletsInGPU.regressionF[quintupletIndex];
@@ -2593,12 +2615,11 @@ namespace SDL {
                                                 lowerModuleIndex4,
                                                 lowerModuleIndex5,
                                                 rPhiChiSquared);
+
       if (not pass)
         return pass;
     }
 
-    float xPix[] = {mdsInGPU.anchorX[pixelInnerMDIndex], mdsInGPU.anchorX[pixelOuterMDIndex]};
-    float yPix[] = {mdsInGPU.anchorY[pixelInnerMDIndex], mdsInGPU.anchorY[pixelOuterMDIndex]};
     rPhiChiSquaredInwards =
         computePT5RPhiChiSquaredInwards(modulesInGPU, T5CenterX, T5CenterY, quintupletRadius, xPix, yPix);
 
@@ -2623,29 +2644,71 @@ namespace SDL {
 
   template <typename TAcc>
   ALPAKA_FN_ACC ALPAKA_FN_INLINE float computePT5RZChiSquared(TAcc const& acc,
-                                                              struct SDL::modules& modulesInGPU,
-                                                              uint16_t* lowerModuleIndices,
-                                                              float* rtPix,
-                                                              float* zPix,
-                                                              float* rts,
-                                                              float* zs) {
-    //use the two anchor hits of the pixel segment to compute the slope
-    //then compute the pseudo chi squared of the five outer hits
-
-    float slope = (zPix[1] - zPix[0]) / (rtPix[1] - rtPix[0]);
+                                                              struct SDL::modules const& modulesInGPU,
+                                                              const uint16_t* lowerModuleIndices,
+                                                              const float* rtPix,
+                                                              const float* xPix,
+                                                              const float* yPix,
+                                                              const float* zPix,
+                                                              const float* rts,
+                                                              const float* zs,
+                                                              float pixelSegmentPt,
+                                                              float pixelSegmentPx,
+                                                              float pixelSegmentPy,
+                                                              float pixelSegmentPz,
+                                                              int pixelSegmentCharge) {
     float residual = 0;
     float error = 0;
-    //hardcoded array indices!!!
     float RMSE = 0;
+
+    // the pixel positions are in unit of cm, and need to be divided by 100 to be in consistent with unit mm.
+    float Px = pixelSegmentPx, Py = pixelSegmentPy, Pz = pixelSegmentPz;
+    int charge = pixelSegmentCharge;
+    float x1 = xPix[1] / 100;
+    float y1 = yPix[1] / 100;
+    float z1 = zPix[1] / 100;
+    float r1 = rtPix[1] / 100;
+
+    float a = -100 / SDL::kR1GeVf * charge;
+
     for (size_t i = 0; i < 5; i++) {
-      uint16_t& lowerModuleIndex = lowerModuleIndices[i];
+      float zsi = zs[i] / 100;
+      float rtsi = rts[i] / 100;
+      uint16_t lowerModuleIndex = lowerModuleIndices[i];
       const int moduleType = modulesInGPU.moduleType[lowerModuleIndex];
       const int moduleSide = modulesInGPU.sides[lowerModuleIndex];
       const int moduleSubdet = modulesInGPU.subdets[lowerModuleIndex];
 
-      residual = (moduleSubdet == SDL::Barrel) ? (zs[i] - zPix[0]) - slope * (rts[i] - rtPix[0])
-                                               : (rts[i] - rtPix[0]) - (zs[i] - zPix[0]) / slope;
-      const float& drdz = modulesInGPU.drdzs[lowerModuleIndex];
+      // calculation is detailed documented here https://indico.cern.ch/event/1185895/contributions/4982756/attachments/2526561/4345805/helix%20pT3%20summarize.pdf
+      float diffr, diffz;
+      float p = alpaka::math::sqrt(acc, Px * Px + Py * Py + Pz * Pz);
+
+      float rou = a / p;
+      if (moduleSubdet == SDL::Endcap) {
+        float s = (zsi - z1) * p / Pz;
+        float x = x1 + Px / a * alpaka::math::sin(acc, rou * s) - Py / a * (1 - alpaka::math::cos(acc, rou * s));
+        float y = y1 + Py / a * alpaka::math::sin(acc, rou * s) + Px / a * (1 - alpaka::math::cos(acc, rou * s));
+        diffr = alpaka::math::abs(acc, rtsi - alpaka::math::sqrt(acc, x * x + y * y)) * 100;
+      }
+
+      if (moduleSubdet == SDL::Barrel) {
+        float paraA = r1 * r1 + 2 * (Px * Px + Py * Py) / (a * a) + 2 * (y1 * Px - x1 * Py) / a - rtsi * rtsi;
+        float paraB = 2 * (x1 * Px + y1 * Py) / a;
+        float paraC = 2 * (y1 * Px - x1 * Py) / a + 2 * (Px * Px + Py * Py) / (a * a);
+        float A = paraB * paraB + paraC * paraC;
+        float B = 2 * paraA * paraB;
+        float C = paraA * paraA - paraC * paraC;
+        float sol1 = (-B + alpaka::math::sqrt(acc, B * B - 4 * A * C)) / (2 * A);
+        float sol2 = (-B - alpaka::math::sqrt(acc, B * B - 4 * A * C)) / (2 * A);
+        float solz1 = alpaka::math::asin(acc, sol1) / rou * Pz / p + z1;
+        float solz2 = alpaka::math::asin(acc, sol2) / rou * Pz / p + z1;
+        float diffz1 = alpaka::math::abs(acc, solz1 - zsi) * 100;
+        float diffz2 = alpaka::math::abs(acc, solz2 - zsi) * 100;
+        diffz = alpaka::math::min(acc, diffz1, diffz2);
+      }
+
+      residual = moduleSubdet == SDL::Barrel ? diffz : diffr;
+
       //PS Modules
       if (moduleType == 0) {
         error = 0.15f;
@@ -2656,12 +2719,14 @@ namespace SDL {
 
       //special dispensation to tilted PS modules!
       if (moduleType == 0 and moduleSubdet == SDL::Barrel and moduleSide != Center) {
-        error /= alpaka::math::sqrt(acc, 1.f + drdz * drdz);
+        float drdz = modulesInGPU.drdzs[lowerModuleIndex];
+        error /= alpaka::math::sqrt(acc, 1 + drdz * drdz);
       }
       RMSE += (residual * residual) / (error * error);
     }
 
-    RMSE = alpaka::math::sqrt(acc, 0.2f * RMSE);
+    RMSE = alpaka::math::sqrt(acc, 0.2f * RMSE);  //Divided by the degree of freedom 5.
+
     return RMSE;
   };
 

RecoTracker/LSTCore/standalone/code/core/write_sdl_ntuple.cc

@@ -404,6 +404,7 @@ void setPixelQuintupletOutputBranches(SDL::Event<SDL::Acc>* event) {
     ana.tx->pushbackToBranch<float>("pT5_phi", phi);
     ana.tx->pushbackToBranch<int>("pT5_layer_binary", layer_binary);
     ana.tx->pushbackToBranch<int>("pT5_moduleType_binary", moduleType_binary);
+    ana.tx->pushbackToBranch<float>("pT5_rzChiSquared", pixelQuintupletsInGPU.rzChiSquared[pT5]);
 
     pT5_matched_simIdx.push_back(simidx);