|
21 | 21 | #include "TPCBase/ParameterGas.h" |
22 | 22 | #include "TPCBase/Sector.h" |
23 | 23 | #include "SpacePoints/TrackResiduals.h" |
| 24 | +#include "SpacePoints/SpacePointsCalibConfParam.h" |
24 | 25 | #include "DataFormatsTPC/Defs.h" |
25 | 26 | #include "TPCFastTransform.h" |
26 | 27 | #include "ChebyshevFit1D.h" |
@@ -513,21 +514,83 @@ std::unique_ptr<o2::gpu::TPCFastSpaceChargeCorrection> TPCFastSpaceChargeCorrect |
513 | 514 | knotsDouble[0].reserve(nY2Xbins); |
514 | 515 | knotsDouble[1].reserve(nZ2Xbins); |
515 | 516 |
|
516 | | - // to get enouth measurements, make a spline knot at every second bin. Boundary bins are always included. |
517 | | - |
518 | | - for (int i = 0, j = nY2Xbins - 1; i <= j; i += 2, j -= 2) { |
519 | | - knotsDouble[0].push_back(trackResiduals.getY2X(0, i)); |
520 | | - if (j >= i + 1) { |
521 | | - knotsDouble[0].push_back(trackResiduals.getY2X(0, j)); |
522 | | - } |
| 517 | + // the spline knots, in the coordinate of the knot calculation below (y/x of row 0, z/x): the bin centres and widths of the voxels |
| 518 | + std::vector<double> binCentres[2]; |
| 519 | + std::vector<double> binWidths[2]; |
| 520 | + for (int i = 0; i < nY2Xbins; i++) { |
| 521 | + binCentres[0].push_back(trackResiduals.getY2X(0, i)); |
| 522 | + binWidths[0].push_back(1. / trackResiduals.getDY2XI(0, i)); |
523 | 523 | } |
| 524 | + for (int i = 0; i < nZ2Xbins; i++) { |
| 525 | + binCentres[1].push_back(trackResiduals.getZ2X(i)); |
| 526 | + binWidths[1].push_back(trackResiduals.getDZ2X(i)); |
| 527 | + } |
| 528 | + // explicit knots (setY2XKnots / setZ2XKnots) are given in the coordinate of the voxel binning: y/x scaled to -1..1, z/x scaled to 0..1 |
| 529 | + const double knotScale[2] = {trackResiduals.getMaxY2X(0), o2::tpc::SpacePointsCalibConfParam::Instance().maxZ2X}; |
| 530 | + const std::vector<double>* explicitKnots[2] = {&mY2XKnots, &mZ2XKnots}; |
524 | 531 |
|
525 | | - for (int i = 0, j = nZ2Xbins - 1; i <= j; i += 2, j -= 2) { |
526 | | - knotsDouble[1].push_back(trackResiduals.getZ2X(i)); |
527 | | - if (j >= i + 1) { |
528 | | - knotsDouble[1].push_back(trackResiduals.getZ2X(j)); |
| 532 | + for (int dim = 0; dim < 2; dim++) { |
| 533 | + const auto& centres = binCentres[dim]; |
| 534 | + auto& knots = knotsDouble[dim]; |
| 535 | + if (explicitKnots[dim]->empty()) { |
| 536 | + // default: to get enough measurements, a spline knot at every second bin, counted from both ends. Boundary bins are always included. |
| 537 | + for (int i = 0, j = int(centres.size()) - 1; i <= j; i += 2, j -= 2) { |
| 538 | + knots.push_back(centres[i]); |
| 539 | + if (j >= i + 1) { |
| 540 | + knots.push_back(centres[j]); |
| 541 | + } |
| 542 | + } |
| 543 | + continue; |
| 544 | + } |
| 545 | + // explicit knots: the spline of a row spans its first to its last bin centre, so those are knots and the others must lie between them. |
| 546 | + // A knot closer than half the narrowest voxel bin to its neighbouring knot is dropped: the knots are put on an integer grid with 1/10 of |
| 547 | + // the smallest knot distance as step, and the spline stores a table with one entry per grid step, so nearly coinciding knots would make |
| 548 | + // the correction object arbitrarily large |
| 549 | + const double first = centres.front(); |
| 550 | + const double last = centres.back(); |
| 551 | + const double eps = 1.e-6 * std::abs(last - first); |
| 552 | + const double minKnotDistance = 0.5 * *std::min_element(binWidths[dim].begin(), binWidths[dim].end()); |
| 553 | + std::vector<double> innerKnots; |
| 554 | + int nOutside = 0; |
| 555 | + for (const double k : *explicitKnots[dim]) { |
| 556 | + const double u = k * knotScale[dim]; |
| 557 | + if (u < first - eps || u > last + eps) { |
| 558 | + ++nOutside; |
| 559 | + } else { |
| 560 | + innerKnots.push_back(u); |
| 561 | + } |
| 562 | + } |
| 563 | + std::sort(innerKnots.begin(), innerKnots.end()); |
| 564 | + int nTooClose = 0; |
| 565 | + knots.push_back(first); |
| 566 | + for (const double u : innerKnots) { |
| 567 | + if (std::abs(u - knots.back()) <= eps || std::abs(u - last) <= eps) { |
| 568 | + continue; // a repeated knot, or the first / the last bin centre given explicitly |
| 569 | + } |
| 570 | + if (u - knots.back() < minKnotDistance || last - u < minKnotDistance) { |
| 571 | + ++nTooClose; |
| 572 | + continue; |
| 573 | + } |
| 574 | + knots.push_back(u); |
| 575 | + } |
| 576 | + knots.push_back(last); |
| 577 | + if (nOutside > 0) { |
| 578 | + LOG(warning) << "fast space charge correction helper: " << nOutside << (dim == 0 ? " y/x" : " z/x") << " knots outside the first and the last bin centre dropped"; |
| 579 | + } |
| 580 | + if (nTooClose > 0) { |
| 581 | + LOG(warning) << "fast space charge correction helper: " << nTooClose << (dim == 0 ? " y/x" : " z/x") << " knots closer than half the narrowest voxel bin (" |
| 582 | + << minKnotDistance / knotScale[dim] << ") to the neighbouring knot dropped"; |
| 583 | + } |
| 584 | + // a knot interval without a voxel centre (endpoints included) leaves the spline there without data |
| 585 | + for (size_t i = 1; i < knots.size(); i++) { |
| 586 | + const bool covered = std::any_of(centres.begin(), centres.end(), [&](double c) { return c >= knots[i - 1] - eps && c <= knots[i] + eps; }); |
| 587 | + if (!covered) { |
| 588 | + LOG(warning) << "fast space charge correction helper: no " << (dim == 0 ? "y/x" : "z/x") << " voxel centre between the knots " << knots[i - 1] / knotScale[dim] << " and " << knots[i] / knotScale[dim]; |
| 589 | + } |
529 | 590 | } |
530 | 591 | } |
| 592 | + LOG(info) << "fast space charge correction helper: " << knotsDouble[0].size() << " y/x knots for " << nY2Xbins << " bins (" << (mY2XKnots.empty() ? "default" : "explicit") << "), " |
| 593 | + << knotsDouble[1].size() << " z/x knots for " << nZ2Xbins << " bins (" << (mZ2XKnots.empty() ? "default" : "explicit") << ")"; |
531 | 594 |
|
532 | 595 | std::vector<int> knotsInt[2]; |
533 | 596 |
|
|
0 commit comments