From d520c207f27ee5bf501638f5e662993846842138 Mon Sep 17 00:00:00 2001 From: YonsiG Date: Tue, 11 Jun 2024 12:17:41 -0700 Subject: [PATCH 1/4] change pT5 rz chi2 cut definition into helix approximation --- RecoTracker/LSTCore/src/alpaka/PixelTriplet.h | 202 ++++++++++++------ .../standalone/code/core/write_sdl_ntuple.cc | 1 + 2 files changed, 133 insertions(+), 70 deletions(-) diff --git a/RecoTracker/LSTCore/src/alpaka/PixelTriplet.h b/RecoTracker/LSTCore/src/alpaka/PixelTriplet.h index bd048f9c819a2..c01b893d1c77e 100644 --- a/RecoTracker/LSTCore/src/alpaka/PixelTriplet.h +++ b/RecoTracker/LSTCore/src/alpaka/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,9 +2556,34 @@ 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; + + //get the appropriate radii and centers + centerX = segmentsInGPU.circleCenterX[pixelSegmentArrayIndex]; + centerY = segmentsInGPU.circleCenterY[pixelSegmentArrayIndex]; + pixelRadius = segmentsInGPU.circleRadius[pixelSegmentArrayIndex]; - if (/*pixelRadius*/ 0 < 5.0f * kR1GeVf) { // FIXME: pixelRadius is not defined yet + 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,11 +2607,6 @@ namespace SDL { mdsInGPU.anchorY[fourthMDIndex], mdsInGPU.anchorY[fifthMDIndex]}; - //get the appropriate radii and centers - centerX = segmentsInGPU.circleCenterX[pixelSegmentArrayIndex]; - centerY = segmentsInGPU.circleCenterY[pixelSegmentArrayIndex]; - pixelRadius = segmentsInGPU.circleRadius[pixelSegmentArrayIndex]; - float T5CenterX = quintupletsInGPU.regressionG[quintupletIndex]; float T5CenterY = quintupletsInGPU.regressionF[quintupletIndex]; quintupletRadius = quintupletsInGPU.regressionRadius[quintupletIndex]; @@ -2603,12 +2622,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 +2651,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; @@ -2666,12 +2726,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/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); From 97e3fb05ebbf697820290b62ea4648e4f6024975 Mon Sep 17 00:00:00 2001 From: YonsiG Date: Wed, 12 Jun 2024 09:23:34 -0700 Subject: [PATCH 2/4] add links and more clear explanations --- RecoTracker/LSTCore/src/alpaka/PixelTriplet.h | 25 +++++++++++-------- 1 file changed, 14 insertions(+), 11 deletions(-) diff --git a/RecoTracker/LSTCore/src/alpaka/PixelTriplet.h b/RecoTracker/LSTCore/src/alpaka/PixelTriplet.h index c01b893d1c77e..162d13ba7f599 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 = -1000/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,6 +2062,8 @@ 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) { // cat 10 return rzChiSquared < 14.031f; @@ -2540,10 +2541,10 @@ 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], @@ -2564,9 +2565,7 @@ namespace SDL { rzChiSquared = 0; - //get the appropriate radii and centers - centerX = segmentsInGPU.circleCenterX[pixelSegmentArrayIndex]; - centerY = segmentsInGPU.circleCenterY[pixelSegmentArrayIndex]; + //get the appropriate centers pixelRadius = segmentsInGPU.circleRadius[pixelSegmentArrayIndex]; if (pixelRadius < 5.0f * kR1GeVf) { //only apply r-z chi2 cuts for <5GeV tracks @@ -2607,6 +2606,10 @@ namespace SDL { mdsInGPU.anchorY[fourthMDIndex], mdsInGPU.anchorY[fifthMDIndex]}; + //get the appropriate centers + centerX = segmentsInGPU.circleCenterX[pixelSegmentArrayIndex]; + centerY = segmentsInGPU.circleCenterY[pixelSegmentArrayIndex]; + float T5CenterX = quintupletsInGPU.regressionG[quintupletIndex]; float T5CenterY = quintupletsInGPU.regressionF[quintupletIndex]; quintupletRadius = quintupletsInGPU.regressionRadius[quintupletIndex]; @@ -2668,6 +2671,7 @@ namespace SDL { float error = 0; 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; @@ -2675,8 +2679,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 = -1000/SDL::kR1GeVf * charge; for (size_t i = 0; i < 5; i++) { float zsi = zs[i] / 100; From 1ccae7ab7a877921a794f33a0d67e4cdbb57b446 Mon Sep 17 00:00:00 2001 From: YonsiG Date: Fri, 14 Jun 2024 14:55:03 -0700 Subject: [PATCH 3/4] bug fix for constant unit --- RecoTracker/LSTCore/src/alpaka/PixelTriplet.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/RecoTracker/LSTCore/src/alpaka/PixelTriplet.h b/RecoTracker/LSTCore/src/alpaka/PixelTriplet.h index 162d13ba7f599..682926f1d8c23 100644 --- a/RecoTracker/LSTCore/src/alpaka/PixelTriplet.h +++ b/RecoTracker/LSTCore/src/alpaka/PixelTriplet.h @@ -760,7 +760,7 @@ namespace SDL { float z1 = zPix[1] / 100; float r1 = rtPix[1] / 100; - float a = -1000/SDL::kR1GeVf * charge; + float a = -100/SDL::kR1GeVf * charge; for (size_t i = 0; i < 3; i++) { float zsi = zs[i] / 100; @@ -2679,7 +2679,7 @@ namespace SDL { float z1 = zPix[1] / 100; float r1 = rtPix[1] / 100; - float a = -1000/SDL::kR1GeVf * charge; + float a = -100/SDL::kR1GeVf * charge; for (size_t i = 0; i < 5; i++) { float zsi = zs[i] / 100; From 5b2fe32cbe94eb106da77b80ce694e8da15ac77c Mon Sep 17 00:00:00 2001 From: YonsiG Date: Mon, 17 Jun 2024 05:20:17 -0700 Subject: [PATCH 4/4] linter checks --- RecoTracker/LSTCore/src/alpaka/PixelTriplet.h | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/RecoTracker/LSTCore/src/alpaka/PixelTriplet.h b/RecoTracker/LSTCore/src/alpaka/PixelTriplet.h index 682926f1d8c23..5a8fa5ed366e9 100644 --- a/RecoTracker/LSTCore/src/alpaka/PixelTriplet.h +++ b/RecoTracker/LSTCore/src/alpaka/PixelTriplet.h @@ -760,7 +760,7 @@ namespace SDL { float z1 = zPix[1] / 100; float r1 = rtPix[1] / 100; - float a = -100/SDL::kR1GeVf * charge; + float a = -100 / SDL::kR1GeVf * charge; for (size_t i = 0; i < 3; i++) { float zsi = zs[i] / 100; @@ -799,7 +799,7 @@ namespace SDL { diffz = alpaka::math::min(acc, diffz1, diffz2); residual = diffz; } - + //PS Modules if (moduleType == 0) { error = 0.15f; @@ -2679,7 +2679,7 @@ namespace SDL { float z1 = zPix[1] / 100; float r1 = rtPix[1] / 100; - float a = -100/SDL::kR1GeVf * charge; + float a = -100 / SDL::kR1GeVf * charge; for (size_t i = 0; i < 5; i++) { float zsi = zs[i] / 100; @@ -2735,7 +2735,7 @@ namespace SDL { RMSE += (residual * residual) / (error * error); } - RMSE = alpaka::math::sqrt(acc, 0.2f * RMSE); //the constant doesn't really matter.... + RMSE = alpaka::math::sqrt(acc, 0.2f * RMSE); //Divided by the degree of freedom 5. return RMSE; };