diff --git a/RecoTracker/LSTCore/src/alpaka/PixelTriplet.h b/RecoTracker/LSTCore/src/alpaka/PixelTriplet.h index bd048f9c819a2..5a8fa5ed366e9 100644 --- a/RecoTracker/LSTCore/src/alpaka/PixelTriplet.h +++ b/RecoTracker/LSTCore/src/alpaka/PixelTriplet.h @@ -760,8 +760,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; @@ -781,6 +780,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) { @@ -797,10 +797,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; @@ -2063,63 +2062,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; @@ -2544,8 +2541,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], @@ -2557,9 +2557,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, @@ -2583,10 +2606,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]; @@ -2603,12 +2625,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); @@ -2633,29 +2654,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; + + // 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; @@ -2666,12 +2729,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; }; diff --git a/RecoTracker/LSTCore/standalone/code/core/write_sdl_ntuple.cc b/RecoTracker/LSTCore/standalone/code/core/write_sdl_ntuple.cc index 66b397434855f..e6eb181cd7dbd 100644 --- a/RecoTracker/LSTCore/standalone/code/core/write_sdl_ntuple.cc +++ b/RecoTracker/LSTCore/standalone/code/core/write_sdl_ntuple.cc @@ -324,6 +324,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);