Skip to content
This repository was archived by the owner on Dec 9, 2024. It is now read-only.

New pT5 r-z chi2 cut definition study based on helix approximation #410

Closed
wants to merge 2 commits into from
Closed

New pT5 r-z chi2 cut definition study based on helix approximation #410

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
dd2a75e
add pT5 rz-chi2 cut threshold helix
YonsiG Jun 11, 2024
2200a3c
linter check
YonsiG Jun 11, 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
189 changes: 125 additions & 64 deletions SDL/PixelTriplet.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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],
Expand All @@ -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
Expand All @@ -2565,7 +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, 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,
Expand Down Expand Up @@ -2604,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);

Expand All @@ -2634,29 +2651,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;

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;
Expand All @@ -2667,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;
};

Expand Down
1 change: 1 addition & 0 deletions code/core/write_sdl_ntuple.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -424,6 +424,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);

Expand Down