From dd2a75e7ba1190d56822016042d8ea74954abae3 Mon Sep 17 00:00:00 2001 From: YonsiG Date: Tue, 11 Jun 2024 07:48:41 -0700 Subject: [PATCH 1/2] add pT5 rz-chi2 cut threshold helix --- SDL/PixelTriplet.h | 177 ++++++++++++++++++++++------------ code/core/write_sdl_ntuple.cc | 1 + 2 files changed, 114 insertions(+), 64 deletions(-) diff --git a/SDL/PixelTriplet.h b/SDL/PixelTriplet.h index 99162649..0ba90adc 100644 --- a/SDL/PixelTriplet.h +++ b/SDL/PixelTriplet.h @@ -2064,62 +2064,58 @@ namespace SDL { modulesInGPU.moduleType[lowerModuleIndex5] == SDL::TwoS); 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; @@ -2546,6 +2542,9 @@ namespace SDL { 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 zs[5] = {mdsInGPU.anchorZ[firstMDIndex], mdsInGPU.anchorZ[secondMDIndex], mdsInGPU.anchorZ[thirdMDIndex], @@ -2557,6 +2556,12 @@ namespace SDL { mdsInGPU.anchorRt[fourthMDIndex], mdsInGPU.anchorRt[fifthMDIndex]}; + 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; //get the appropriate radii and centers @@ -2565,7 +2570,8 @@ namespace SDL { pixelRadius = segmentsInGPU.circleRadius[pixelSegmentArrayIndex]; if (pixelRadius < 5.0f * kR1GeVf) { //only apply r-z chi2 cuts for <5GeV tracks - rzChiSquared = computePT5RZChiSquared(acc, modulesInGPU, lowerModuleIndices, rtPix, zPix, rts, zs); + rzChiSquared = computePT5RZChiSquared(acc, modulesInGPU, lowerModuleIndices, rtPix, xPix, yPix, zPix, rts, zs, + pixelSegmentPt, pixelSegmentPx, pixelSegmentPy, pixelSegmentPz, pixelSegmentCharge); pass = pass and passPT5RZChiSquaredCuts(modulesInGPU, lowerModuleIndex1, lowerModuleIndex2, @@ -2604,12 +2610,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); @@ -2634,29 +2639,71 @@ namespace SDL { template 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; + + 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 Bz = SDL::magnetic_field; + float a = -0.299792 * Bz * 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; @@ -2667,12 +2714,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); //the constant doesn't really matter.... + return RMSE; }; diff --git a/code/core/write_sdl_ntuple.cc b/code/core/write_sdl_ntuple.cc index dad0c4bf..a4599157 100644 --- a/code/core/write_sdl_ntuple.cc +++ b/code/core/write_sdl_ntuple.cc @@ -424,6 +424,7 @@ void setPixelQuintupletOutputBranches(SDL::Event* event) ana.tx->pushbackToBranch("pT5_phi", phi); ana.tx->pushbackToBranch("pT5_layer_binary", layer_binary); ana.tx->pushbackToBranch("pT5_moduleType_binary", moduleType_binary); + ana.tx->pushbackToBranch("pT5_rzChiSquared", pixelQuintupletsInGPU.rzChiSquared[pT5]); pT5_matched_simIdx.push_back(simidx); From 2200a3cf34e585f2af97da6ce2cfc7c112328374 Mon Sep 17 00:00:00 2001 From: YonsiG Date: Tue, 11 Jun 2024 07:53:22 -0700 Subject: [PATCH 2/2] linter check --- SDL/PixelTriplet.h | 56 ++++++++++++++++++++++++++++------------------ 1 file changed, 34 insertions(+), 22 deletions(-) diff --git a/SDL/PixelTriplet.h b/SDL/PixelTriplet.h index 0ba90adc..c01b893d 100644 --- a/SDL/PixelTriplet.h +++ b/SDL/PixelTriplet.h @@ -2064,57 +2064,57 @@ namespace SDL { modulesInGPU.moduleType[lowerModuleIndex5] == SDL::TwoS); if (layer1 == 1 and layer2 == 2 and layer3 == 3) { - if (layer4 == 12 and layer5 == 13) { // cat 10 + if (layer4 == 12 and layer5 == 13) { // cat 10 return rzChiSquared < 14.031f; - } else if (layer4 == 4 and layer5 == 12) { // cat 12 + } else if (layer4 == 4 and layer5 == 12) { // cat 12 return rzChiSquared < 8.760f; - } else if (layer4 == 4 and layer5 == 5) { // cat 11 + } else if (layer4 == 4 and layer5 == 5) { // cat 11 return rzChiSquared < 3.607f; - } else if (layer4 == 7 and layer5 == 13) { // cat 9 + } else if (layer4 == 7 and layer5 == 13) { // cat 9 return rzChiSquared < 16.620; - } else if (layer4 == 7 and layer5 == 8) { // cat 8 + } 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) { // cat 7 + if (layer4 == 13 and layer5 == 14) { // cat 7 return rzChiSquared < 8.950f; - } else if (layer4 == 8 and layer5 == 14) { // cat 6 + } else if (layer4 == 8 and layer5 == 14) { // cat 6 return rzChiSquared < 14.837f; - } else if (layer4 == 8 and layer5 == 9) { // cat 5 + } 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) { // cat 3 + if (layer4 == 9 and layer5 == 10) { // cat 3 return rzChiSquared < 15.093f; - } else if (layer4 == 9 and layer5 == 15) { // cat 4 + } 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) { // cat 20 + if (layer4 == 12 and layer5 == 13) { // cat 20 return rzChiSquared < 12.868f; - } else if (layer4 == 5 and layer5 == 12) { // cat 19 + } else if (layer4 == 5 and layer5 == 12) { // cat 19 return rzChiSquared < 6.128f; - } else if (layer4 == 5 and layer5 == 6) { // cat 18 + } 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) { // cat 17 + if (layer4 == 13 and layer5 == 14) { // cat 17 return rzChiSquared < 19.446f; - } else if (layer4 == 8 and layer5 == 14) { // cat 16 + } 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) { // cat 15 + if (layer4 == 14 and layer5 == 15) { // cat 15 return rzChiSquared < 14.71f; - } else if (layer4 == 9 and layer5 == 15) { // cat 14 + } else if (layer4 == 9 and layer5 == 15) { // cat 14 return rzChiSquared < 18.213f; } } else if (layer1 == 7 and layer2 == 8 and layer3 == 9) { - if (layer4 == 10 and layer5 == 11) { // cat 0 + if (layer4 == 10 and layer5 == 11) { // cat 0 return rzChiSquared < 10.016f; - } else if (layer4 == 10 and layer5 == 16) { // cat 1 + } else if (layer4 == 10 and layer5 == 16) { // cat 1 return rzChiSquared < 87.671f; - } else if (layer4 == 15 and layer5 == 16) { // cat 2 + } else if (layer4 == 15 and layer5 == 16) { // cat 2 return rzChiSquared < 5.844f; } } @@ -2570,8 +2570,20 @@ namespace SDL { 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); + rzChiSquared = computePT5RZChiSquared(acc, + modulesInGPU, + lowerModuleIndices, + rtPix, + xPix, + yPix, + zPix, + rts, + zs, + pixelSegmentPt, + pixelSegmentPx, + pixelSegmentPy, + pixelSegmentPz, + pixelSegmentCharge); pass = pass and passPT5RZChiSquaredCuts(modulesInGPU, lowerModuleIndex1, lowerModuleIndex2,