Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

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

Merged
merged 4 commits into from
Jun 19, 2024
Merged

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

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
d520c20
change pT5 rz chi2 cut definition into helix approximation
YonsiG Jun 11, 2024
97e3fb0
add links and more clear explanations
YonsiG Jun 12, 2024
1ccae7a
bug fix for constant unit
YonsiG Jun 14, 2024
5b2fe32
linter checks
YonsiG Jun 17, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
209 changes: 137 additions & 72 deletions RecoTracker/LSTCore/src/alpaka/PixelTriplet.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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) {
Expand All @@ -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;
Expand Down Expand Up @@ -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
YonsiG marked this conversation as resolved.
Show resolved Hide resolved
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;
Expand Down Expand Up @@ -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]};
YonsiG marked this conversation as resolved.
Show resolved Hide resolved
float zPix[2] = {mdsInGPU.anchorZ[pixelInnerMDIndex], mdsInGPU.anchorZ[pixelOuterMDIndex]};

float zs[5] = {mdsInGPU.anchorZ[firstMDIndex],
mdsInGPU.anchorZ[secondMDIndex],
mdsInGPU.anchorZ[thirdMDIndex],
Expand All @@ -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,
Expand All @@ -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];
Expand All @@ -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);

Expand All @@ -2633,29 +2654,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;
YonsiG marked this conversation as resolved.
Show resolved Hide resolved
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;
YonsiG marked this conversation as resolved.
Show resolved Hide resolved

//PS Modules
if (moduleType == 0) {
error = 0.15f;
Expand All @@ -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;
};

Expand Down
1 change: 1 addition & 0 deletions RecoTracker/LSTCore/standalone/code/core/write_sdl_ntuple.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -324,6 +324,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]);
VourMa marked this conversation as resolved.
Show resolved Hide resolved

pT5_matched_simIdx.push_back(simidx);

Expand Down