Skip to content

Commit 66a90f9

Browse files
cbmswdavidrohr
cbmsw
authored and
davidrohr
committed
TPC Splines: better smoothing between the voxels
1 parent df5e021 commit 66a90f9

File tree

5 files changed

+94
-60
lines changed

5 files changed

+94
-60
lines changed

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

Lines changed: 65 additions & 39 deletions
Original file line numberDiff line numberDiff line change
@@ -148,9 +148,13 @@ void TPCFastSpaceChargeCorrectionHelper::fillSpaceChargeCorrectionFromMap(TPCFas
148148
if (!processingInverseCorrection) {
149149
info.resetMaxValues();
150150
}
151+
info.updateMaxValues(1., 1., 1.);
152+
info.updateMaxValues(-1., -1., -1.);
153+
151154
if (nDataPoints >= 4) {
152155
std::vector<double> pointGU(nDataPoints);
153156
std::vector<double> pointGV(nDataPoints);
157+
std::vector<double> pointWeight(nDataPoints);
154158
std::vector<double> pointCorr(3 * nDataPoints); // 3 dimensions
155159
for (int i = 0; i < nDataPoints; ++i) {
156160
o2::gpu::TPCFastSpaceChargeCorrectionMap::CorrectionPoint p = data[i];
@@ -161,15 +165,14 @@ void TPCFastSpaceChargeCorrectionHelper::fillSpaceChargeCorrectionFromMap(TPCFas
161165
}
162166
pointGU[i] = gu;
163167
pointGV[i] = gv;
168+
pointWeight[i] = p.mWeight;
164169
pointCorr[3 * i + 0] = p.mDx;
165170
pointCorr[3 * i + 1] = p.mDy;
166171
pointCorr[3 * i + 2] = p.mDz;
167-
if (!processingInverseCorrection) {
168-
info.updateMaxValues(20. * p.mDx, 20. * p.mDy, 20. * p.mDz);
169-
}
172+
info.updateMaxValues(5. * p.mDx, 5. * p.mDy, 5. * p.mDz);
170173
}
171-
helper.approximateDataPoints(spline, splineParameters.data(), 0., spline.getGridX1().getUmax(), 0., spline.getGridX2().getUmax(), &pointGU[0],
172-
&pointGV[0], &pointCorr[0], nDataPoints);
174+
helper.approximateDataPoints(spline, splineParameters.data(), 0., spline.getGridX1().getUmax(), 0., spline.getGridX2().getUmax(), pointGU.data(),
175+
pointGV.data(), pointCorr.data(), pointWeight.data(), nDataPoints);
173176
} else {
174177
for (int i = 0; i < spline.getNumberOfParameters(); i++) {
175178
splineParameters[i] = 0.f;
@@ -301,7 +304,7 @@ std::unique_ptr<TPCFastSpaceChargeCorrection> TPCFastSpaceChargeCorrectionHelper
301304
double dx, dy, dz;
302305
correctionLocal(iSector, iRow, y, z, dx, dy, dz);
303306
mCorrectionMap.addCorrectionPoint(iSector, iRow,
304-
y, z, dx, dy, dz);
307+
y, z, dx, dy, dz, 1.);
305308
}
306309
}
307310
} // row
@@ -593,14 +596,6 @@ std::unique_ptr<o2::gpu::TPCFastSpaceChargeCorrection> TPCFastSpaceChargeCorrect
593596
data.mCy *= -1.;
594597
data.mCz *= -1.;
595598
}
596-
if (data.mNentries > 0) {
597-
if (iSector < geo.getNumberOfSectorsA() && data.mZ < 0) {
598-
LOG(error) << errMsg << "fitted Z coordinate " << data.mZ << " is negative for sector " << iSector;
599-
}
600-
if (iSector >= geo.getNumberOfSectorsA() && data.mZ > 0) {
601-
LOG(error) << errMsg << "fitted Z coordinate " << data.mZ << " is positive for sector " << iSector;
602-
}
603-
}
604599
}
605600
};
606601
processor.Process(myThread);
@@ -624,6 +619,9 @@ std::unique_ptr<o2::gpu::TPCFastSpaceChargeCorrection> TPCFastSpaceChargeCorrect
624619
}
625620
}
626621

622+
double maxError[3] = {0., 0., 0.};
623+
int nErrors = 0;
624+
627625
for (int iSector = 0; iSector < nSectors; iSector++) {
628626

629627
// now process the data row-by-row
@@ -682,9 +680,21 @@ std::unique_ptr<o2::gpu::TPCFastSpaceChargeCorrection> TPCFastSpaceChargeCorrect
682680
}
683681

684682
if (!msg.str().empty()) {
685-
LOG(warning) << directionName << " correction: fitted voxel position is outside the voxel: "
686-
<< " sector " << iSector << " row " << iRow << " bin y " << iy << " bin z " << iz
687-
<< msg.str();
683+
bool isMaxErrorExceeded = (fabs(data.mX - x) / dx > maxError[0]) ||
684+
(fabs(data.mY - vox.mY) / vox.mDy > maxError[1]) ||
685+
(fabs(data.mZ - vox.mZ) / vox.mDz > maxError[2]);
686+
static std::mutex mutex;
687+
mutex.lock();
688+
nErrors++;
689+
if (nErrors < 20 || isMaxErrorExceeded) {
690+
LOG(warning) << directionName << " correction: error N " << nErrors << "fitted voxel position is outside the voxel: "
691+
<< " sector " << iSector << " row " << iRow << " bin y " << iy << " bin z " << iz
692+
<< msg.str();
693+
maxError[0] = GPUCommonMath::Max(maxError[0], fabs(data.mX - x) / dx);
694+
maxError[1] = GPUCommonMath::Max(maxError[1], fabs(data.mY - vox.mY) / vox.mDy);
695+
maxError[2] = GPUCommonMath::Max(maxError[2], fabs(data.mZ - vox.mZ) / vox.mDz);
696+
}
697+
mutex.unlock();
688698
}
689699

690700
} else { // no data, take voxel center position
@@ -773,17 +783,27 @@ std::unique_ptr<o2::gpu::TPCFastSpaceChargeCorrection> TPCFastSpaceChargeCorrect
773783
auto& info = correction.getSectorRowInfo(iSector, iRow);
774784
const auto& spline = correction.getSpline(iSector, iRow);
775785

776-
auto addEdge = [&](int iy1, int iz1, int iy2, int iz2, int nSteps) {
786+
auto addVoxel = [&](int iy, int iz, double weight) {
787+
auto& vox = vRowVoxels[iy * nZ2Xbins + iz];
788+
if (vox.mSmoothingStep > 2) {
789+
LOG(fatal) << "empty voxel is not repared: y " << iy << " z " << iz;
790+
}
791+
auto& data = vSectorData[iSector * nRows + iRow][iy * nZ2Xbins + iz];
792+
map.addCorrectionPoint(iSector, iRow, data.mY, data.mZ, data.mCx, data.mCy, data.mCz, weight);
793+
};
794+
795+
auto addEdge = [&](int iy1, int iz1, int iy2, int iz2, int nPoints) {
796+
// add n points on the edge between two voxels excluding the voxel points
797+
if (nPoints < 1)
798+
return;
799+
if (iy1 < 0 || iy1 >= nY2Xbins || iz1 < 0 || iz1 >= nZ2Xbins)
800+
return;
801+
if (iy2 < 0 || iy2 >= nY2Xbins || iz2 < 0 || iz2 >= nZ2Xbins)
802+
return;
777803
auto& data1 = vSectorData[iSector * nRows + iRow][iy1 * nZ2Xbins + iz1];
778804
auto& vox1 = vRowVoxels[iy1 * nZ2Xbins + iz1];
779805
auto& data2 = vSectorData[iSector * nRows + iRow][iy2 * nZ2Xbins + iz2];
780806
auto& vox2 = vRowVoxels[iy2 * nZ2Xbins + iz2];
781-
if (vox1.mSmoothingStep > 2) {
782-
LOG(fatal) << "empty voxel is not repared: y " << iy1 << " z " << iz1;
783-
}
784-
if (vox2.mSmoothingStep > 2) {
785-
LOG(fatal) << "empty voxel is not repared: y " << iy2 << " z " << iz2;
786-
}
787807
double y1 = data1.mY;
788808
double z1 = data1.mZ;
789809
double cx1 = data1.mCx;
@@ -795,32 +815,36 @@ std::unique_ptr<o2::gpu::TPCFastSpaceChargeCorrection> TPCFastSpaceChargeCorrect
795815
double cy2 = data2.mCy;
796816
double cz2 = data2.mCz;
797817

798-
for (int is = 0; is < nSteps; is++) {
799-
double s2 = is / (double)nSteps;
818+
for (int is = 1; is <= nPoints; is++) {
819+
double s2 = is / (double)(nPoints + 1);
800820
double s1 = 1. - s2;
801821
double y = s1 * y1 + s2 * y2;
802822
double z = s1 * z1 + s2 * z2;
803823
double cx = s1 * cx1 + s2 * cx2;
804824
double cy = s1 * cy1 + s2 * cy2;
805825
double cz = s1 * cz1 + s2 * cz2;
806-
map.addCorrectionPoint(iSector, iRow, y, z, cx, cy, cz);
826+
map.addCorrectionPoint(iSector, iRow, y, z, cx, cy, cz, 1.);
807827
}
808828
};
809829

830+
// original measurements weighted by 8 at each voxel and 8 additional artificial measurements around each voxel
831+
//
832+
// (y+1, z) 8 1 1 8 (y+1, z+1)
833+
// 1 1 1 1 1
834+
// 1 1 1 1 1
835+
// (y,z) 8 1 1 8 1
836+
// 1 1 1 1 1
837+
810838
for (int iy = 0; iy < nY2Xbins; iy++) {
811-
for (int iz = 0; iz < nZ2Xbins - 1; iz++) {
812-
addEdge(iy, iz, iy, iz + 1, 3);
839+
for (int iz = 0; iz < nZ2Xbins; iz++) {
840+
addVoxel(iy, iz, 8);
841+
addEdge(iy, iz, iy, iz + 1, 2);
842+
addEdge(iy, iz, iy + 1, iz, 2);
843+
addEdge(iy, iz, iy + 1, iz + 1, 2);
844+
addEdge(iy + 1, iz, iy, iz + 1, 2);
813845
}
814-
addEdge(iy, nZ2Xbins - 1, iy, nZ2Xbins - 1, 1);
815846
}
816847

817-
for (int iz = 0; iz < nZ2Xbins; iz++) {
818-
for (int iy = 0; iy < nY2Xbins - 1; iy++) {
819-
addEdge(iy, iz, iy + 1, iz, 3);
820-
}
821-
addEdge(nY2Xbins - 1, iz, nY2Xbins - 1, iz, 1);
822-
} // iy
823-
824848
} // iRow
825849
}; // myThread
826850

@@ -939,10 +963,11 @@ void TPCFastSpaceChargeCorrectionHelper::initInverse(std::vector<o2::gpu::TPCFas
939963
}
940964
}
941965

942-
std::vector<double> dataPointGridU, dataPointGridV, dataPointF;
966+
std::vector<double> dataPointGridU, dataPointGridV, dataPointF, dataPointWeight;
943967
dataPointGridU.reserve(gridU.size() * gridV.size());
944968
dataPointGridV.reserve(gridU.size() * gridV.size());
945969
dataPointF.reserve(3 * gridU.size() * gridV.size());
970+
dataPointWeight.reserve(gridU.size() * gridV.size());
946971

947972
for (int iu = 0; iu < gridU.size(); iu++) {
948973
for (int iv = 0; iv < gridV.size(); iv++) {
@@ -963,6 +988,7 @@ void TPCFastSpaceChargeCorrectionHelper::initInverse(std::vector<o2::gpu::TPCFas
963988
dataPointF.push_back(scale * dx);
964989
dataPointF.push_back(scale * dy);
965990
dataPointF.push_back(scale * dz);
991+
dataPointWeight.push_back(1.);
966992
}
967993
}
968994

@@ -972,7 +998,7 @@ void TPCFastSpaceChargeCorrectionHelper::initInverse(std::vector<o2::gpu::TPCFas
972998
helper.approximateDataPoints(spline, splineParameters.data(), 0., spline.getGridX1().getUmax(),
973999
0., spline.getGridX2().getUmax(),
9741000
dataPointGridU.data(), dataPointGridV.data(),
975-
dataPointF.data(), nDataPoints);
1001+
dataPointF.data(), dataPointWeight.data(), nDataPoints);
9761002

9771003
float* splineX = correction.getSplineDataInvX(sector, row);
9781004
float* splineUV = correction.getSplineDataInvYZ(sector, row);

‎GPU/TPCFastTransformation/Spline2DHelper.cxx‎

Lines changed: 10 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -241,6 +241,7 @@ void Spline2DHelper<DataT>::approximateFunctionViaDataPoints(
241241
mFdimensions = spline.getYdimensions();
242242
std::vector<double> dataPointX1(getNumberOfDataPoints());
243243
std::vector<double> dataPointX2(getNumberOfDataPoints());
244+
std::vector<double> dataPointWeight(getNumberOfDataPoints(), 1.);
244245
std::vector<double> dataPointF(getNumberOfDataPoints() * mFdimensions);
245246

246247
double scaleX1 = (x1Max - x1Min) / ((double)mHelperU1.getSpline().getUmax());
@@ -256,7 +257,8 @@ void Spline2DHelper<DataT>::approximateFunctionViaDataPoints(
256257
F(x1, x2, &dataPointF[ind * mFdimensions]);
257258
}
258259
}
259-
approximateDataPoints(spline, spline.getParameters(), x1Min, x1Max, x2Min, x2Max, &dataPointX1[0], &dataPointX2[0], &dataPointF[0], getNumberOfDataPoints());
260+
approximateDataPoints(spline, spline.getParameters(), x1Min, x1Max, x2Min, x2Max, dataPointX1.data(), dataPointX2.data(), dataPointF.data(),
261+
dataPointWeight.data(), getNumberOfDataPoints());
260262
}
261263

262264
template <typename DataT>
@@ -326,7 +328,7 @@ template <typename DataT>
326328
void Spline2DHelper<DataT>::approximateDataPoints(
327329
Spline2DContainer<DataT>& spline, DataT* splineParameters, double x1Min, double x1Max, double x2Min, double x2Max,
328330
const double dataPointX1[], const double dataPointX2[], const double dataPointF[/*getNumberOfDataPoints() x nFdim*/],
329-
int32_t nDataPoints)
331+
const double dataPointWeight[], int32_t nDataPoints)
330332
{
331333
/// Create best-fit spline parameters for a given input function F
332334

@@ -343,6 +345,10 @@ void Spline2DHelper<DataT>::approximateDataPoints(
343345
for (int32_t iPoint = 0; iPoint < nDataPoints; ++iPoint) {
344346
double u = fGridU.convXtoU(dataPointX1[iPoint]);
345347
double v = fGridV.convXtoU(dataPointX2[iPoint]);
348+
double weight = dataPointWeight[iPoint];
349+
if (!(weight > 0.)) {
350+
continue;
351+
}
346352
int32_t iu = fGridU.getLeftKnotIndexForU(u);
347353
int32_t iv = fGridV.getLeftKnotIndexForU(v);
348354
double c[16];
@@ -353,14 +359,14 @@ void Spline2DHelper<DataT>::approximateDataPoints(
353359

354360
for (int32_t i = 0; i < 16; i++) {
355361
for (int32_t j = i; j < 16; j++) {
356-
solver.A(ind[i], ind[j]) += c[i] * c[j];
362+
solver.A(ind[i], ind[j]) += weight * c[i] * c[j];
357363
}
358364
}
359365

360366
for (int32_t iDim = 0; iDim < nFdim; iDim++) {
361367
double f = (double)dataPointF[iPoint * nFdim + iDim];
362368
for (int32_t i = 0; i < 16; i++) {
363-
solver.B(ind[i], iDim) += f * c[i];
369+
solver.B(ind[i], iDim) += weight * f * c[i];
364370
}
365371
}
366372
} // data points

‎GPU/TPCFastTransformation/Spline2DHelper.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -74,7 +74,7 @@ class Spline2DHelper
7474
void approximateDataPoints(
7575
Spline2DContainer<DataT>& spline, DataT* splineParameters, double x1Min, double x1Max, double x2Min, double x2Max,
7676
const double dataPointX1[/*nDataPoints*/], const double dataPointX2[/*nDataPoints*/],
77-
const double dataPointF[/*nDataPoints x spline.getYdimensions*/], int32_t nDataPoints);
77+
const double dataPointF[/*nDataPoints x spline.getYdimensions*/], const double dataPointWeight[/*nDataPoints*/], int32_t nDataPoints);
7878

7979
/// _______________ Interface for a step-wise construction of the best-fit spline ________________________
8080

‎GPU/TPCFastTransformation/TPCFastSpaceChargeCorrectionMap.h‎

Lines changed: 5 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -42,8 +42,9 @@ class TPCFastSpaceChargeCorrectionMap
4242
/// \brief The struct contains necessary info for TPC padrow
4343
///
4444
struct CorrectionPoint {
45-
double mY, mZ; // not-distorted local coordinates
46-
double mDx, mDy, mDz; // corrections to the local coordinates
45+
double mY{0.}, mZ{0.}; // not-distorted local coordinates
46+
double mDx{0.}, mDy{0.}, mDz{0.}; // corrections to the local coordinates
47+
double mWeight{0.}; // weight of the point
4748
};
4849

4950
/// _____________ Constructors / destructors __________________________
@@ -72,11 +73,11 @@ class TPCFastSpaceChargeCorrectionMap
7273
/// Starts the construction procedure, reserves temporary memory
7374
void addCorrectionPoint(int32_t iSector, int32_t iRow,
7475
double y, double z,
75-
double dx, double dy, double dz)
76+
double dx, double dy, double dz, double weight)
7677
{
7778
int32_t ind = mNrows * iSector + iRow;
7879
fDataPoints.at(ind).push_back(CorrectionPoint{y, z,
79-
dx, dy, dz});
80+
dx, dy, dz, weight});
8081
}
8182

8283
const std::vector<CorrectionPoint>& getPoints(int32_t iSector, int32_t iRow) const

‎GPU/TPCFastTransformation/macro/TPCFastTransformInit.C‎

Lines changed: 13 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -319,6 +319,8 @@ void TPCFastTransformInit(const char* fileName = "debugVoxRes.root", const char*
319319

320320
for (int direction = 0; direction < 2; direction++) { // 0 - normal, 1 - inverse
321321

322+
std::string directionName = (direction == 0) ? "direct" : "inverse";
323+
322324
TTree* currentTree = (direction == 0) ? voxResTree : voxResTreeInverse;
323325
if (!currentTree) {
324326
std::cout << "tree voxResTree does not exist!" << std::endl;
@@ -342,7 +344,7 @@ void TPCFastTransformInit(const char* fileName = "debugVoxRes.root", const char*
342344
double sumDiff[3] = {0., 0., 0.};
343345
int64_t nDiff = 0;
344346

345-
std::cout << "fill debug ntuples at voxels ..." << std::endl;
347+
LOG(info) << directionName << " correction: fill debug ntuples at voxels ...";
346348

347349
for (int32_t iVox = 0; iVox < currentTree->GetEntriesFast(); iVox++) {
348350

@@ -457,12 +459,11 @@ void TPCFastTransformInit(const char* fileName = "debugVoxRes.root", const char*
457459
}
458460
}
459461

460-
std::cout
461-
<< "fill debug ntuples everywhere .." << std::endl;
462+
LOG(info) << directionName << " correction: fill debug ntuples everywhere ..";
462463

463464
for (int32_t iSector = 0; iSector < geo.getNumberOfSectors(); iSector++) {
464465
// for (int32_t iSector = 0; iSector < 1; iSector++) {
465-
std::cout << "debug ntules for sector " << iSector << std::endl;
466+
LOG(info) << directionName << " correction: fill debug ntuples everywhere in sector " << iSector;
466467

467468
int mirrorSector = (iSector >= geo.getNumberOfSectorsA()) ? iSector - geo.getNumberOfSectorsA() : iSector + geo.getNumberOfSectorsA();
468469

@@ -592,17 +593,17 @@ void TPCFastTransformInit(const char* fileName = "debugVoxRes.root", const char*
592593
for (int32_t i = 0; i < 3; i++) {
593594
sumDiff[i] = sqrt(sumDiff[i]) / nDiff;
594595
}
596+
LOG(info) << directionName << " correction: max and mean differences between spline and voxel corrections:";
597+
LOG(info) << "Max difference in x : " << maxDiff[0] << " at Sector "
598+
<< maxDiffSector[0] << " row " << maxDiffRow[0];
595599

596-
std::cout << "Max difference in x : " << maxDiff[0] << " at Sector "
597-
<< maxDiffSector[0] << " row " << maxDiffRow[0] << std::endl;
598-
599-
std::cout << "Max difference in y : " << maxDiff[1] << " at Sector "
600-
<< maxDiffSector[1] << " row " << maxDiffRow[1] << std::endl;
600+
LOG(info) << "Max difference in y : " << maxDiff[1] << " at Sector "
601+
<< maxDiffSector[1] << " row " << maxDiffRow[1];
601602

602-
std::cout << "Max difference in z : " << maxDiff[2] << " at Sector "
603-
<< maxDiffSector[2] << " row " << maxDiffRow[2] << std::endl;
603+
LOG(info) << "Max difference in z : " << maxDiff[2] << " at Sector "
604+
<< maxDiffSector[2] << " row " << maxDiffRow[2];
604605

605-
std::cout << "Mean difference in x,y,z : " << sumDiff[0] << " " << sumDiff[1]
606+
LOG(info) << "Mean difference in x,y,z : " << sumDiff[0] << " " << sumDiff[1]
606607
<< " " << sumDiff[2] << std::endl;
607608
} // direction
608609

0 commit comments

Comments
 (0)