diff --git a/GPU/TPCFastTransformation/TPCFastSpaceChargeCorrection.cxx b/GPU/TPCFastTransformation/TPCFastSpaceChargeCorrection.cxx index a8f46373d2ff5..9230034354b62 100644 --- a/GPU/TPCFastTransformation/TPCFastSpaceChargeCorrection.cxx +++ b/GPU/TPCFastTransformation/TPCFastSpaceChargeCorrection.cxx @@ -147,9 +147,8 @@ void TPCFastSpaceChargeCorrection::setActualBufferAddress(char* actualFlatBuffer } size_t bufferSize = scBufferOffset + scBufferSize; for (int32_t is = 0; is < 3; is++) { - size_t correctionDataOffset = alignSize(bufferSize, SplineType::getParameterAlignmentBytes()); - mCorrectionData[is] = reinterpret_cast(mFlatBufferPtr + correctionDataOffset); - bufferSize = correctionDataOffset + mSectorDataSizeBytes[is] * mGeo.getNumberOfSectors(); + mCorrectionData[is] = reinterpret_cast(mFlatBufferPtr + bufferSize); + bufferSize += mSectorDataSizeBytes[is] * mGeo.getNumberOfSectors(); } return; } @@ -256,7 +255,7 @@ void TPCFastSpaceChargeCorrection::setActualBufferAddress(char* actualFlatBuffer for (int32_t is = 0; is < 3; is++) { size_t oldCorrectionDataOffset = alignSize(oldBufferSize, SplineType::getParameterAlignmentBytes()); - size_t correctionDataOffset = alignSize(bufferSize, SplineType::getParameterAlignmentBytes()); + size_t correctionDataOffset = bufferSize; mCorrectionData[is] = reinterpret_cast(mFlatBufferPtr + correctionDataOffset); memmove(mCorrectionData[is], mFlatBufferPtr + oldCorrectionDataOffset, mSectorDataSizeBytes[is] * mGeo.getNumberOfSectors()); oldBufferSize = oldCorrectionDataOffset + mSectorDataSizeBytes[is] * mGeo.getNumberOfSectors(); @@ -509,16 +508,24 @@ void TPCFastSpaceChargeCorrection::finishConstruction() size_t bufferSize = scBufferOffsets[0] + scBufferSize; size_t correctionDataOffset[3]; for (int32_t is = 0; is < 3; is++) { - correctionDataOffset[is] = alignSize(bufferSize, SplineType::getParameterAlignmentBytes()); + correctionDataOffset[is] = bufferSize; mSectorDataSizeBytes[is] = 0; for (int32_t j = 0; j < mGeo.getNumberOfRows(); j++) { RowInfo& row = getRowInfo(j); - SplineType& spline = mConstructionScenarios[row.splineScenarioID]; - row.dataOffsetBytes[is] = alignSize(mSectorDataSizeBytes[is], SplineType::getParameterAlignmentBytes()); - mSectorDataSizeBytes[is] = row.dataOffsetBytes[is] + spline.getSizeOfParameters(); + row.dataOffsetBytes[is] = mSectorDataSizeBytes[is]; + const SplineType& spline = mConstructionScenarios[row.splineScenarioID]; + if (is == 0) { + const SplineTypeXYZ& splineXYZ = reinterpret_cast(spline); + mSectorDataSizeBytes[is] += splineXYZ.getSizeOfParameters(); + } else if (is == 1) { + const SplineTypeInvX& splineInvX = reinterpret_cast(spline); + mSectorDataSizeBytes[is] += splineInvX.getSizeOfParameters(); + } else if (is == 2) { + const SplineTypeInvYZ& splineInvYZ = reinterpret_cast(spline); + mSectorDataSizeBytes[is] += splineInvYZ.getSizeOfParameters(); + } } - mSectorDataSizeBytes[is] = alignSize(mSectorDataSizeBytes[is], SplineType::getParameterAlignmentBytes()); - bufferSize = correctionDataOffset[is] + mSectorDataSizeBytes[is] * mGeo.getNumberOfSectors(); + bufferSize += mSectorDataSizeBytes[is] * mGeo.getNumberOfSectors(); } FlatObject::finishConstruction(bufferSize); @@ -558,24 +565,23 @@ GPUd() void TPCFastSpaceChargeCorrection::setNoCorrection() } // row for (int32_t sector = 0; sector < mGeo.getNumberOfSectors(); sector++) { - for (int32_t row = 0; row < mGeo.getNumberOfRows(); row++) { - const SplineType& spline = getSplineForRow(row); for (int32_t is = 0; is < 3; is++) { float* data = getCorrectionData(sector, row, is); - int32_t nPar = spline.getNumberOfParameters(); - if (is == 1) { - nPar = nPar / 3; - } - if (is == 2) { - nPar = nPar * 2 / 3; + int32_t nPar = 0; + if (is == 0) { + nPar = getSplineForRow(row).getNumberOfParameters(); + } else if (is == 1) { + nPar = getSplineInvXforRow(row).getNumberOfParameters(); + } else if (is == 2) { + nPar = getSplineInvYZforRow(row).getNumberOfParameters(); } for (int32_t i = 0; i < nPar; i++) { data[i] = 0.f; } } } // row - } // sector + } // sector } void TPCFastSpaceChargeCorrection::constructWithNoCorrection(const TPCFastTransformGeo& geo) diff --git a/GPU/TPCFastTransformation/TPCFastTransformPOD.cxx b/GPU/TPCFastTransformation/TPCFastTransformPOD.cxx index 0709b96dad0a9..f6e61bdbbff70 100644 --- a/GPU/TPCFastTransformation/TPCFastTransformPOD.cxx +++ b/GPU/TPCFastTransformation/TPCFastTransformPOD.cxx @@ -31,7 +31,7 @@ TPCFastTransformPOD* TPCFastTransformPOD::create(aligned_unique_buffer_ptrsetTimeStamp(src.getTimeStamp()); res->setVDrift(src.getVDrift()); @@ -48,7 +48,7 @@ TPCFastTransformPOD* TPCFastTransformPOD::create(aligned_unique_buffer_ptr(buff + nextDynOffs); LOGP(debug, "splinID={} start offset {} -> {}", is, nextDynOffs, (void*)data); // metadata - size_t sectorDataSizeFloats = 0; - size_t sectorDataSizeBytes = 0; - for (int row = 0; row < NROWS; row++) { - podMap.getRowInfo(row).dataOffsetBytes[is] = sectorDataSizeBytes; - const auto& spline = origCorr.getSplineForRow(row); - int nPar = spline.getNumberOfParameters(); - if (is == 1) { - nPar = nPar / 3; - } - if (is == 2) { - nPar = nPar * 2 / 3; - } - sectorDataSizeFloats += nPar; - sectorDataSizeBytes += nPar * sizeof(float); - } + size_t sectorDataSizeBytes = origCorr.mSectorDataSizeBytes[is]; - for (int sector = 0; sector < origCorr.mGeo.getNumberOfSectors(); sector++) { - podMap.mSplineDataOffsets[sector][is] = nextDynOffs; - if (buffSize < nextDynOffs + sectorDataSizeBytes) { - throw std::runtime_error(fmt::format("attempt to copy {} bytes of data for spline{} of sector{} to {}, overflowing the buffer of size {}", sectorDataSizeBytes, is, sector, nextDynOffs, buffSize)); - } - const float* dataOr = origCorr.getCorrectionData(sector, 0, is); - std::memcpy(data, dataOr, sectorDataSizeBytes); - data += sectorDataSizeFloats; - nextDynOffs += sectorDataSizeBytes; + for (int sector = 0; sector < TPCFastTransformGeo::getNumberOfSectors(); sector++) { + podMap.mSplineDataOffsets[sector][is] = nextDynOffs + sectorDataSizeBytes * sector; } + size_t dataSize = TPCFastTransformGeo::getNumberOfSectors() * sectorDataSizeBytes; + if (buffSize < nextDynOffs + dataSize) { + throw std::runtime_error(fmt::format("attempt to copy {} bytes of data for spline{} to {}, overflowing the buffer of size {}", sectorDataSizeBytes, is, nextDynOffs, buffSize)); + } + const char* dataOr = origCorr.mCorrectionData[is]; + std::memcpy(data, dataOr, dataSize); + nextDynOffs += dataSize; } - podMap.mTotalSize = alignOffset(nextDynOffs); + nextDynOffs = alignOffset(nextDynOffs); + podMap.mTotalSize = nextDynOffs; if (buffSize != podMap.mTotalSize) { throw std::runtime_error(fmt::format("Estimated buffer size {} differs from filled one {}", buffSize, podMap.mTotalSize)); } diff --git a/GPU/TPCFastTransformation/macro/TPCFastTransformInit.C b/GPU/TPCFastTransformation/macro/TPCFastTransformInit.C index dd8bb1a789db9..dace124617cac 100644 --- a/GPU/TPCFastTransformation/macro/TPCFastTransformInit.C +++ b/GPU/TPCFastTransformation/macro/TPCFastTransformInit.C @@ -210,6 +210,14 @@ void TPCFastTransformInit(const char* fileName = "debugVoxRes.root", const char* o2::gpu::TPCFastSpaceChargeCorrection& corr = fastTransform->getCorrection(); + aligned_unique_buffer_ptr podBuffer; + TPCFastTransformPOD* corrPODptr = TPCFastTransformPOD::create(podBuffer, *fastTransform); + + if (!corrPODptr) { + throw std::runtime_error("Failed to create TPCFastTransformPOD"); + } + const TPCFastTransformPOD& corrPOD = *corrPODptr; + // a debug file with some NTuples TDirectory* currDir = gDirectory; @@ -302,19 +310,35 @@ void TPCFastTransformInit(const char* fileName = "debugVoxRes.root", const char* const o2::gpu::TPCFastTransformGeo& geo = helper->getGeometry(); + float maxPodDiff[6] = {0., 0., 0., 0., 0., 0.}; + auto getInvCorrections = [&](int iSector, int iRow, float realY, float realZ, float& ix, float& iy, float& iz) { // get the inverse corrections ix, iy, iz at x,y,z ix = corr.getCorrectionXatRealYZ(iSector, iRow, realY, realZ); corr.getCorrectionYZatRealYZ(iSector, iRow, realY, realZ, iy, iz); + + float ixPod = corrPOD.getCorrectionXatRealYZ(iSector, iRow, realY, realZ); + float iyPod, izPod; + corrPOD.getCorrectionYZatRealYZ(iSector, iRow, realY, realZ, iyPod, izPod); + + maxPodDiff[3] = std::max(maxPodDiff[3], fabs(ix - ixPod)); + maxPodDiff[4] = std::max(maxPodDiff[4], fabs(iy - iyPod)); + maxPodDiff[5] = std::max(maxPodDiff[5], fabs(iz - izPod)); }; auto getAllCorrections = [&](int iSector, int iRow, float y, float z, float& cx, float& cy, float& cz, float& ix, float& iy, float& iz) { // get the corrections cx,cy,cz at x,y,z corr.getCorrectionLocal(iSector, iRow, y, z, cx, cy, cz); getInvCorrections(iSector, iRow, y + cy, z + cz, ix, iy, iz); + + float cxPod, cyPod, czPod; + corrPOD.getCorrectionLocal(iSector, iRow, y, z, cxPod, cyPod, czPod); + maxPodDiff[0] = std::max(maxPodDiff[0], fabs(cx - cxPod)); + maxPodDiff[1] = std::max(maxPodDiff[1], fabs(cy - cyPod)); + maxPodDiff[2] = std::max(maxPodDiff[2], fabs(cz - czPod)); }; - for (int direction = 0; direction < 2; direction++) { // 0 - normal, 1 - inverse + for (int direction = 0; direction < 2; direction++) { // 0 - direct, 1 - inverse std::string directionName = (direction == 0) ? "direct" : "inverse"; @@ -594,7 +618,9 @@ void TPCFastTransformInit(const char* fileName = "debugVoxRes.root", const char* for (int32_t i = 0; i < 3; i++) { sumDiff[i] = sqrt(sumDiff[i]) / nDiff; } + LOG(info) << directionName << " correction: max and mean differences between spline and voxel corrections:"; + LOG(info) << "Max difference in x : " << maxDiff[0] << " at Sector " << maxDiffSector[0] << " row " << maxDiffRow[0]; @@ -605,9 +631,22 @@ void TPCFastTransformInit(const char* fileName = "debugVoxRes.root", const char* << maxDiffSector[2] << " row " << maxDiffRow[2]; LOG(info) << "Mean difference in x,y,z : " << sumDiff[0] << " " << sumDiff[1] - << " " << sumDiff[2] << std::endl; + << " " << sumDiff[2]; + + LOG(info) << std::endl; + } // direction + LOG(info) << " max difference between POD and original corrections: "; + LOG(info) << " x " << maxPodDiff[0]; + LOG(info) << " y " << maxPodDiff[1]; + LOG(info) << " z " << maxPodDiff[2]; + LOG(info) << " inverse x " << maxPodDiff[3]; + LOG(info) << " inverse y " << maxPodDiff[4]; + LOG(info) << " inverse z " << maxPodDiff[5]; + + LOG(info) << std::endl; + corr.testInverse(true); debugFile->cd();