Skip to content

Commit cdb2b5f

Browse files
Sergey Gorbunovcbmsw
Sergey Gorbunov
authored and
cbmsw
committed
TPC Splines: smooth to linear edges, crop at grid borders, use mean position of residuals
1 parent 7416b3e commit cdb2b5f

File tree

4 files changed

+172
-183
lines changed

4 files changed

+172
-183
lines changed

‎Detectors/TPC/calibration/src/TPCFastSpaceChargeCorrectionHelper.cxx‎

Lines changed: 127 additions & 175 deletions
Original file line numberDiff line numberDiff line change
@@ -159,8 +159,8 @@ void TPCFastSpaceChargeCorrectionHelper::fillSpaceChargeCorrectionFromMap(TPCFas
159159
pointCorr[3 * i + 0] = dx;
160160
pointCorr[3 * i + 1] = du;
161161
pointCorr[3 * i + 2] = dv;
162-
info.updateMaxValues(2. * dx, 2. * du, 2. * dv);
163-
info.updateMaxValuesInv(-2. * dx, -2. * du, -2. * dv);
162+
info.updateMaxValues(20. * dx, 20. * du, 20. * dv);
163+
info.updateMaxValuesInv(-20. * dx, -20. * du, -20. * dv);
164164
}
165165
helper.approximateDataPoints(spline, splineParameters, 0., spline.getGridX1().getUmax(), 0., spline.getGridX2().getUmax(), &pointSU[0],
166166
&pointSV[0], &pointCorr[0], nDataPoints);
@@ -411,95 +411,69 @@ std::unique_ptr<o2::gpu::TPCFastSpaceChargeCorrection> TPCFastSpaceChargeCorrect
411411
int nY2Xbins = trackResiduals.getNY2XBins();
412412
int nZ2Xbins = trackResiduals.getNZ2XBins();
413413

414-
double marginY2X = trackResiduals.getY2X(0, 2) - trackResiduals.getY2X(0, 0);
415-
double marginZ2X = trackResiduals.getZ2X(1) - trackResiduals.getZ2X(0);
414+
std::vector<double> uvBinsDouble[2];
416415

417-
std::vector<int> yBinsInt;
418-
{
419-
std::vector<double> yBins;
420-
yBins.reserve(nY2Xbins + 2);
421-
yBins.push_back(trackResiduals.getY2X(0, 0) - marginY2X);
422-
for (int i = 0, j = nY2Xbins - 1; i <= j; i += 2, j -= 2) {
423-
if (i == j) {
424-
yBins.push_back(trackResiduals.getY2X(0, i));
425-
} else if (i + 1 == j) {
426-
yBins.push_back(trackResiduals.getY2X(0, i));
427-
} else {
428-
yBins.push_back(trackResiduals.getY2X(0, i));
429-
yBins.push_back(trackResiduals.getY2X(0, j));
430-
}
416+
uvBinsDouble[0].reserve(nY2Xbins);
417+
uvBinsDouble[1].reserve(nZ2Xbins);
418+
419+
for (int i = 0, j = nY2Xbins - 1; i <= j; i += 2, j -= 2) {
420+
uvBinsDouble[0].push_back(trackResiduals.getY2X(0, i));
421+
if (j >= i + 1) {
422+
uvBinsDouble[0].push_back(trackResiduals.getY2X(0, j));
431423
}
432-
yBins.push_back(trackResiduals.getY2X(0, nY2Xbins - 1) + marginY2X);
424+
}
433425

434-
std::sort(yBins.begin(), yBins.end());
435-
double dy = yBins[1] - yBins[0];
436-
for (int i = 1; i < yBins.size(); i++) {
437-
if (yBins[i] - yBins[i - 1] < dy) {
438-
dy = yBins[i] - yBins[i - 1];
426+
for (int i = 0, j = nZ2Xbins - 1; i <= j; i += 2, j -= 2) {
427+
uvBinsDouble[1].push_back(-trackResiduals.getZ2X(i));
428+
if (j >= i + 1) {
429+
uvBinsDouble[1].push_back(-trackResiduals.getZ2X(j));
430+
}
431+
}
432+
433+
std::vector<int> uvBinsInt[2];
434+
435+
for (int iuv = 0; iuv < 2; iuv++) {
436+
auto& bins = uvBinsDouble[iuv];
437+
std::sort(bins.begin(), bins.end());
438+
439+
auto& binsInt = uvBinsInt[iuv];
440+
binsInt.reserve(bins.size());
441+
442+
double dy = bins[1] - bins[0];
443+
for (int i = 2; i < bins.size(); i++) {
444+
double dd = bins[i] - bins[i - 1];
445+
if (dd < dy) {
446+
dy = dd;
439447
}
440448
}
441-
yBinsInt.reserve(yBins.size());
442449
// spline knots must be positioned on the grid with integer internal coordinate
443450
// take the knot position accuracy of 0.1*dy
444451
dy = dy / 10.;
445-
double y0 = yBins[0];
446-
double y1 = yBins[yBins.size() - 1];
447-
for (auto& y : yBins) {
452+
double y0 = bins[0];
453+
double y1 = bins[bins.size() - 1];
454+
for (auto& y : bins) {
448455
y -= y0;
449456
int iy = int(y / dy + 0.5);
450-
yBinsInt.push_back(iy);
457+
binsInt.push_back(iy);
451458
double yold = y / (y1 - y0) * 2 - 1.;
452459
y = iy * dy;
453460
y = y / (y1 - y0) * 2 - 1.;
454-
LOG(info) << "convert y bin: " << yold << " -> " << y << " -> " << iy;
455-
}
456-
}
457-
458-
std::vector<int> zBinsInt;
459-
{
460-
std::vector<double> zBins;
461-
zBins.reserve(nZ2Xbins + 2);
462-
zBins.push_back(-(trackResiduals.getZ2X(0) - marginZ2X));
463-
for (int i = 0; i < nZ2Xbins; i += 2) {
464-
zBins.push_back(-trackResiduals.getZ2X(i));
465-
}
466-
zBins.push_back(-(trackResiduals.getZ2X(nZ2Xbins - 1) + 2. * marginZ2X));
467-
468-
std::sort(zBins.begin(), zBins.end());
469-
double dz = zBins[1] - zBins[0];
470-
for (int i = 1; i < zBins.size(); i++) {
471-
if (zBins[i] - zBins[i - 1] < dz) {
472-
dz = zBins[i] - zBins[i - 1];
461+
if (iuv == 0) {
462+
LOG(info) << "convert y bin: " << yold << " -> " << y << " -> " << iy;
463+
} else {
464+
LOG(info) << "convert z bin: " << yold << " -> " << y << " -> " << iy;
473465
}
474466
}
475-
zBinsInt.reserve(zBins.size());
476-
// spline knots must be positioned on the grid with an integer internal coordinate
477-
// lets copy the knot positions with the accuracy of 0.01*dz
478-
dz = dz / 10.;
479-
double z0 = zBins[0];
480-
double z1 = zBins[zBins.size() - 1];
481-
for (auto& z : zBins) {
482-
z -= z0;
483-
int iz = int(z / dz + 0.5);
484-
zBinsInt.push_back(iz);
485-
double zold = z / (z1 - z0);
486-
z = iz * dz;
487-
z = z / (z1 - z0);
488-
LOG(info) << "convert z bin: " << zold << " -> " << z << " -> " << iz;
489-
}
490-
}
491467

492-
if (yBinsInt.size() < 2) {
493-
yBinsInt.clear();
494-
yBinsInt.push_back(0);
495-
yBinsInt.push_back(1);
468+
if (binsInt.size() < 2) {
469+
binsInt.clear();
470+
binsInt.push_back(0);
471+
binsInt.push_back(1);
472+
}
496473
}
497474

498-
if (zBinsInt.size() < 2) {
499-
zBinsInt.clear();
500-
zBinsInt.push_back(0);
501-
zBinsInt.push_back(1);
502-
}
475+
auto& yBinsInt = uvBinsInt[0];
476+
auto& zBinsInt = uvBinsInt[1];
503477

504478
int nKnotsY = yBinsInt.size();
505479
int nKnotsZ = zBinsInt.size();
@@ -534,10 +508,10 @@ std::unique_ptr<o2::gpu::TPCFastSpaceChargeCorrection> TPCFastSpaceChargeCorrect
534508
const auto& rowInfo = geo.getRowInfo(iRow);
535509
auto& info = correction.getSliceRowInfo(iRoc, iRow);
536510
const auto& spline = correction.getSpline(iRoc, iRow);
537-
double yMin = rowInfo.x * (trackResiduals.getY2X(iRow, 0) - marginY2X);
538-
double yMax = rowInfo.x * (trackResiduals.getY2X(iRow, trackResiduals.getNY2XBins() - 1) + marginY2X);
539-
double zMin = rowInfo.x * (trackResiduals.getZ2X(0) - marginZ2X);
540-
double zMax = rowInfo.x * (trackResiduals.getZ2X(trackResiduals.getNZ2XBins() - 1) + 2. * marginZ2X);
511+
double yMin = rowInfo.x * trackResiduals.getY2X(iRow, 0);
512+
double yMax = rowInfo.x * trackResiduals.getY2X(iRow, trackResiduals.getNY2XBins() - 1);
513+
double zMin = rowInfo.x * trackResiduals.getZ2X(0);
514+
double zMax = rowInfo.x * trackResiduals.getZ2X(trackResiduals.getNZ2XBins() - 1);
541515
double uMin = yMin;
542516
double uMax = yMax;
543517
double vMin = geo.getTPCzLength(iRoc) - zMax;
@@ -563,6 +537,7 @@ std::unique_ptr<o2::gpu::TPCFastSpaceChargeCorrection> TPCFastSpaceChargeCorrect
563537

564538
struct VoxelData {
565539
int mNentries{0}; // number of entries
540+
float mX, mY, mZ; // mean position in the local coordinates
566541
float mCx, mCy, mCz; // corrections to the local coordinates
567542
};
568543

@@ -589,16 +564,19 @@ std::unique_ptr<o2::gpu::TPCFastSpaceChargeCorrection> TPCFastSpaceChargeCorrect
589564
}
590565
int iy = v->bvox[o2::tpc::TrackResiduals::VoxF]; // bin number in y/x 0..14
591566
int iz = v->bvox[o2::tpc::TrackResiduals::VoxZ]; // bin number in z/x 0..4
592-
auto& vox = vRocData[iRoc * nRows + iRow][iy * nZ2Xbins + iz];
593-
vox.mNentries = (int)v->stat[o2::tpc::TrackResiduals::VoxV];
594-
vox.mCx = useSmoothed ? v->DS[o2::tpc::TrackResiduals::ResX] : v->D[o2::tpc::TrackResiduals::ResX];
595-
vox.mCy = useSmoothed ? v->DS[o2::tpc::TrackResiduals::ResY] : v->D[o2::tpc::TrackResiduals::ResY];
596-
vox.mCz = useSmoothed ? v->DS[o2::tpc::TrackResiduals::ResZ] : v->D[o2::tpc::TrackResiduals::ResZ];
597-
if (0 && vox.mNentries < 1) {
598-
vox.mCx = 0.;
599-
vox.mCy = 0.;
600-
vox.mCz = 0.;
601-
vox.mNentries = 1;
567+
auto& data = vRocData[iRoc * nRows + iRow][iy * nZ2Xbins + iz];
568+
data.mNentries = (int)v->stat[o2::tpc::TrackResiduals::VoxV];
569+
data.mX = v->stat[o2::tpc::TrackResiduals::VoxX];
570+
data.mY = v->stat[o2::tpc::TrackResiduals::VoxF];
571+
data.mZ = v->stat[o2::tpc::TrackResiduals::VoxZ];
572+
data.mCx = useSmoothed ? v->DS[o2::tpc::TrackResiduals::ResX] : v->D[o2::tpc::TrackResiduals::ResX];
573+
data.mCy = useSmoothed ? v->DS[o2::tpc::TrackResiduals::ResY] : v->D[o2::tpc::TrackResiduals::ResY];
574+
data.mCz = useSmoothed ? v->DS[o2::tpc::TrackResiduals::ResZ] : v->D[o2::tpc::TrackResiduals::ResZ];
575+
if (0 && data.mNentries < 1) {
576+
data.mCx = 0.;
577+
data.mCy = 0.;
578+
data.mCz = 0.;
579+
data.mNentries = 1;
602580
}
603581
}
604582
};
@@ -642,10 +620,27 @@ std::unique_ptr<o2::gpu::TPCFastSpaceChargeCorrection> TPCFastSpaceChargeCorrect
642620
if (iRoc >= geo.getNumberOfSlicesA()) {
643621
vox.mZ = -vox.mZ;
644622
}
623+
data.mY *= x;
624+
data.mZ *= x;
625+
/*
626+
if ( fabs(x - data.mX) > 0.01 || fabs(vox.mY - data.mY) > 5. || fabs(vox.mZ - data.mZ) > 5.) {
627+
std::cout
628+
<< " roc " << iRoc << " row " << iRow
629+
<< " voxel x " << x << " y " << vox.mY << " z " << vox.mZ
630+
<< " data x " << data.mX << " y " << data.mY << " z " << data.mZ
631+
<< std::endl;
632+
}
633+
*/
634+
if (1) { // always use voxel center instead of the mean position
635+
data.mY = vox.mY;
636+
data.mZ = vox.mZ;
637+
}
645638
if (data.mNentries < 1) { // no data
646639
data.mCx = 0.;
647640
data.mCy = 0.;
648641
data.mCz = 0.;
642+
data.mY = vox.mY;
643+
data.mZ = vox.mZ;
649644
vox.mSmoothingStep = 100;
650645
} else { // voxel contains data
651646
if (invertSigns) {
@@ -726,102 +721,59 @@ std::unique_ptr<o2::gpu::TPCFastSpaceChargeCorrection> TPCFastSpaceChargeCorrect
726721

727722
// feed the row data to the helper
728723

729-
double yMin = 0., yMax = 0., zMin = 0.;
730-
731724
auto& info = correction.getSliceRowInfo(iRoc, iRow);
732725
const auto& spline = correction.getSpline(iRoc, iRow);
733726

734-
{
735-
float u0, u1, v0, v1;
736-
correction.convGridToUV(iRoc, iRow, 0., 0., u0, v0);
737-
correction.convGridToUV(iRoc, iRow,
738-
spline.getGridX1().getUmax(), spline.getGridX2().getUmax(), u1, v1);
739-
float y0, y1, z0, z1;
740-
geo.convUVtoLocal(iRoc, u0, v0, y0, z0);
741-
geo.convUVtoLocal(iRoc, u1, v1, y1, z1);
742-
if (iRoc < geo.getNumberOfSlicesA()) {
743-
yMin = y0;
744-
yMax = y1;
745-
} else {
746-
yMin = y1;
747-
yMax = y0;
727+
auto addEdge = [&](int iy1, int iz1, int iy2, int iz2, int nSteps) {
728+
auto& data1 = vRocData[iRoc * nRows + iRow][iy1 * nZ2Xbins + iz1];
729+
auto& vox1 = vRowVoxels[iy1 * nZ2Xbins + iz1];
730+
auto& data2 = vRocData[iRoc * nRows + iRow][iy2 * nZ2Xbins + iz2];
731+
auto& vox2 = vRowVoxels[iy2 * nZ2Xbins + iz2];
732+
if (vox1.mSmoothingStep > 2) {
733+
LOG(fatal) << "empty voxel is not repared: y " << iy1 << " z " << iz1;
748734
}
749-
zMin = z1;
750-
}
751-
752-
double zEdge = 0.;
753-
if (iRoc < geo.getNumberOfSlicesA()) {
754-
zEdge = geo.getTPCzLengthA();
755-
} else {
756-
zEdge = -geo.getTPCzLengthC();
757-
}
735+
if (vox2.mSmoothingStep > 2) {
736+
LOG(fatal) << "empty voxel is not repared: y " << iy2 << " z " << iz2;
737+
}
738+
double y1 = vox1.mY;
739+
double z1 = vox1.mZ;
740+
double cx1 = data1.mCx;
741+
double cy1 = data1.mCy;
742+
double cz1 = data1.mCz;
743+
double y2 = vox2.mY;
744+
double z2 = vox2.mZ;
745+
double cx2 = data2.mCx;
746+
double cy2 = data2.mCy;
747+
double cz2 = data2.mCz;
748+
749+
for (int is = 0; is < nSteps; is++) {
750+
double s2 = is / (double)nSteps;
751+
double s1 = 1. - s2;
752+
double y = s1 * y1 + s2 * y2;
753+
double z = s1 * z1 + s2 * z2;
754+
double cx = s1 * cx1 + s2 * cx2;
755+
double cy = s1 * cy1 + s2 * cy2;
756+
double cz = s1 * cz1 + s2 * cz2;
757+
map.addCorrectionPoint(iRoc, iRow, y, z, cx, cy, cz);
758+
}
759+
};
758760

759761
for (int iy = 0; iy < nY2Xbins; iy++) {
760-
for (int iz = 0; iz < nZ2Xbins; iz++) {
761-
auto& data = vRocData[iRoc * nRows + iRow][iy * nZ2Xbins + iz];
762-
auto& vox = vRowVoxels[iy * nZ2Xbins + iz];
763-
if (vox.mSmoothingStep > 2) {
764-
LOG(fatal) << "empty voxel is not repared";
765-
}
766-
767-
double y = vox.mY;
768-
double z = vox.mZ;
769-
double dy = vox.mDy;
770-
double dz = vox.mDz;
771-
double correctionX = data.mCx;
772-
double correctionY = data.mCy;
773-
double correctionZ = data.mCz;
774-
775-
double yStep = dy / 2.;
776-
double zStep = dz / 2.;
777-
778-
double yFirst = y;
779-
double yLast = y;
780-
double zFirst = z;
781-
double zLast = z;
782-
783-
if (iy == 0) { // extend value of the first Y bin to the row edge
784-
yFirst = yMin;
785-
yStep = (yLast - yFirst) / 2.;
786-
}
787-
788-
if (iy == nY2Xbins - 1) { // extend value of the last Y bin to the row edge
789-
yLast = yMax;
790-
yStep = (yLast - yFirst) / 2.;
791-
}
792-
793-
for (double py = yFirst; py <= yLast + yStep / 2.; py += yStep) {
794-
795-
for (double pz = zFirst; pz <= zLast + zStep / 2.; pz += zStep) {
796-
map.addCorrectionPoint(iRoc, iRow, py, pz, correctionX, correctionY,
797-
correctionZ);
798-
}
762+
for (int iz = 0; iz < nZ2Xbins - 1; iz++) {
763+
addEdge(iy, iz, iy, iz + 1, 3);
764+
}
765+
addEdge(iy, nZ2Xbins - 1, iy, nZ2Xbins - 1, 1);
766+
}
799767

800-
if (iz == 0) { // extend value of the first Z bin to Z=0.
801-
int nZsteps = 2;
802-
for (int is = 0; is < nZsteps; is++) {
803-
double pz = z + (zMin - z) * (is + 1.) / nZsteps;
804-
double s = 1.; //(nZsteps - 1. - is) / nZsteps;
805-
map.addCorrectionPoint(iRoc, iRow, py, pz, s * correctionX,
806-
s * correctionY, s * correctionZ);
807-
}
808-
}
768+
for (int iz = 0; iz < nZ2Xbins; iz++) {
769+
for (int iy = 0; iy < nY2Xbins - 1; iy++) {
770+
addEdge(iy, iz, iy + 1, iz, 3);
771+
}
772+
addEdge(nY2Xbins - 1, iz, nY2Xbins - 1, iz, 1);
773+
} // iy
809774

810-
if (iz == nZ2Xbins - 1) {
811-
// extend value of the last Z bin to the readout, linear decrease of all values to 0.
812-
int nZsteps = 2;
813-
for (int is = 0; is < nZsteps; is++) {
814-
double pz = z + (zEdge - z) * (is + 1.) / nZsteps;
815-
double s = (nZsteps - 1. - is) / nZsteps;
816-
map.addCorrectionPoint(iRoc, iRow, py, pz, s * correctionX,
817-
s * correctionY, s * correctionZ);
818-
}
819-
}
820-
}
821-
} // iz
822-
} // iy
823-
} // iRow
824-
}; // myThread
775+
} // iRow
776+
}; // myThread
825777

826778
// run n threads
827779

0 commit comments

Comments
 (0)