From cd9df99e7cc4e8b37babec6a839110a90c15b5ad Mon Sep 17 00:00:00 2001 From: shahoian Date: Sun, 22 Feb 2026 12:43:33 +0100 Subject: [PATCH] Add POD version of TPCFastTransform The TPCFastTransformPOD is a pointerless version of the TPCFastTransform. It can be created from the original TPCFastTransform as e.g. auto lold = o2::gpu::TPCFastTransform::loadFromFile("o2-gpu-TPCFastTransform.root","ccdb_object"); // load original transform std::vector v; // one has to provide a vector (could be a std or pmr), which later can be messaged via DPL auto* pod = o2::gpu::TPCFastTransformPOD::create(v, *lold); // pointer pod is just v.data() cast to TPCFastTransformPOD* // run test: pod->test(*lold); [INFO] (ns per call) original this Nmissmatch [INFO] getCorrection 1.330e+02 1.400e+02 0 [INFO] getCorrectionInvCorrectedX 8.856e+01 8.434e+01 0 [INFO] getCorrectionInvUV 6.266e+01 6.142e+01 0 It can be also created directly from the TPCFastSpaceChargeCorrection as TPCFastSpaceChargeCorrection& oldCorr = lold->getCorrection(); auto* pod = o2::gpu::TPCFastTransformPOD::create(v, oldCorr); but in this case one should afterwards set the vdrift and t0 using provided getters. TPCFastTransformPOD replicates all the methods of the TPCFastTransform (and of the TPCFastSpaceChargeCorrection), including those which allow to query rescaled corrections (by providing refernce maps and scaling coefficients). Since the idea of this class is to create a final correction map as a weighted sum of different contribution and to distribute it to consumer processes via shared memory, also the query methods w/o rescaling are added, they have the suffix _new added. Eventually, the scalable legacy methods can be suppressed and the suffix new can be dropped. --- GPU/TPCFastTransformation/CMakeLists.txt | 1 + .../TPCFastSpaceChargeCorrection.h | 2 + .../TPCFastTransformPOD.cxx | 238 +++++ .../TPCFastTransformPOD.h | 920 ++++++++++++++++++ .../TPCFastTransformationLinkDef_O2.h | 1 + 5 files changed, 1162 insertions(+) create mode 100644 GPU/TPCFastTransformation/TPCFastTransformPOD.cxx create mode 100644 GPU/TPCFastTransformation/TPCFastTransformPOD.h diff --git a/GPU/TPCFastTransformation/CMakeLists.txt b/GPU/TPCFastTransformation/CMakeLists.txt index 182a66fb28296..769e9981102ef 100644 --- a/GPU/TPCFastTransformation/CMakeLists.txt +++ b/GPU/TPCFastTransformation/CMakeLists.txt @@ -26,6 +26,7 @@ set(SRCS TPCFastSpaceChargeCorrectionMap.cxx TPCFastTransform.cxx CorrectionMapsHelper.cxx + TPCFastTransformPOD.cxx ) if(NOT ALIGPU_BUILD_TYPE STREQUAL "Standalone") diff --git a/GPU/TPCFastTransformation/TPCFastSpaceChargeCorrection.h b/GPU/TPCFastTransformation/TPCFastSpaceChargeCorrection.h index b1a3d0c35da7c..1a11c94888067 100644 --- a/GPU/TPCFastTransformation/TPCFastSpaceChargeCorrection.h +++ b/GPU/TPCFastTransformation/TPCFastSpaceChargeCorrection.h @@ -38,6 +38,8 @@ namespace gpu /// class TPCFastSpaceChargeCorrection : public FlatObject { + friend class TPCFastTransformPOD; + public: // obsolete structure, declared here only for backward compatibility struct SliceInfo { diff --git a/GPU/TPCFastTransformation/TPCFastTransformPOD.cxx b/GPU/TPCFastTransformation/TPCFastTransformPOD.cxx new file mode 100644 index 0000000000000..8796e96f228d5 --- /dev/null +++ b/GPU/TPCFastTransformation/TPCFastTransformPOD.cxx @@ -0,0 +1,238 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file TPCFastTransformPOD.cxx +/// \brief Implementation of POD correction map +/// +/// \author ruben.shahoayn@cern.ch + +#include "TPCFastTransformPOD.h" +/// \brief Implementation of POD correction map +/// +/// \author ruben.shahoayn@cern.ch + +#include "TPCFastTransformPOD.h" +#include "GPUDebugStreamer.h" +#if !defined(GPUCA_GPUCODE) +#include +#endif + +namespace o2 +{ +namespace gpu +{ + +#if !defined(GPUCA_GPUCODE) + +size_t TPCFastTransformPOD::estimateSize(const TPCFastSpaceChargeCorrection& origCorr) +{ + // estimate size of own buffer + const size_t selfSizeFix = sizeof(TPCFastTransformPOD); + size_t nextDynOffs = alignOffset(selfSizeFix); + nextDynOffs = alignOffset(nextDynOffs + origCorr.mNumberOfScenarios * sizeof(size_t)); // spline scenarios start here + // space for splines + for (int isc = 0; isc < origCorr.mNumberOfScenarios; isc++) { + const auto& spline = origCorr.mScenarioPtr[isc]; + nextDynOffs = alignOffset(nextDynOffs + sizeof(spline)); + } + // space for splines data + for (int is = 0; is < 3; is++) { + for (int sector = 0; sector < origCorr.mGeo.getNumberOfSectors(); sector++) { + for (int row = 0; row < NROWS; row++) { + const auto& spline = origCorr.getSpline(sector, row); + int nPar = spline.getNumberOfParameters(); + if (is == 1) { + nPar = nPar / 3; + } + if (is == 2) { + nPar = nPar * 2 / 3; + } + nextDynOffs += nPar * sizeof(float); + } + } + } + nextDynOffs = alignOffset(nextDynOffs); + return nextDynOffs; +} + +TPCFastTransformPOD* TPCFastTransformPOD::create(char* buff, size_t buffSize, const TPCFastSpaceChargeCorrection& origCorr) +{ + // instantiate object to already created buffer of the right size + assert(buffSize > sizeof(TPCFastTransformPOD)); + auto& podMap = getNonConst(buff); + podMap.mApplyCorrection = true; // by default always apply corrections + + // copy fixed size data --- start + podMap.mNumberOfScenarios = origCorr.mNumberOfScenarios; + std::memcpy(&podMap.mGeo, &origCorr.mGeo, sizeof(TPCFastTransformGeo)); // copy geometry (fixed size) + for (int sector = 0; sector < TPCFastTransformGeo::getNumberOfSectors(); sector++) { + for (int row = 0; row < NROWS; row++) { + podMap.mSectorRowInfos[NROWS * sector + row] = origCorr.getSectorRowInfo(sector, row); + } + } + podMap.mTimeStamp = origCorr.mTimeStamp; + // + // init data members coming from the TPCFastTrasform + podMap.mVdrift = 0.; + podMap.mT0 = 0.; + // copy fixed size data --- end + + size_t nextDynOffs = alignOffset(sizeof(TPCFastTransformPOD)); + + // copy sector scenarios + podMap.mOffsScenariosOffsets = nextDynOffs; // spline scenarios offsets start here + LOGP(debug, "Set mOffsScenariosOffsets = {}", podMap.mOffsScenariosOffsets); + nextDynOffs = alignOffset(nextDynOffs + podMap.mNumberOfScenarios * sizeof(size_t)); // spline scenarios start here + + // copy spline objects + size_t* scenOffs = reinterpret_cast(buff + podMap.mOffsScenariosOffsets); + for (int isc = 0; isc < origCorr.mNumberOfScenarios; isc++) { + scenOffs[isc] = nextDynOffs; + const auto& spline = origCorr.mScenarioPtr[isc]; + if (buffSize < nextDynOffs + sizeof(spline)) { + throw std::runtime_error(fmt::format("attempt to copy {} bytes for spline for scenario {} to {}, overflowing the buffer of size {}", sizeof(spline), isc, nextDynOffs + sizeof(spline), buffSize)); + } + std::memcpy(buff + scenOffs[isc], &spline, sizeof(spline)); + nextDynOffs = alignOffset(nextDynOffs + sizeof(spline)); + LOGP(debug, "Copy {} bytes for spline scenario {} (ptr:{}) to offsset {}", sizeof(spline), isc, (void*)&spline, scenOffs[isc]); + } + + // copy splines data + for (int is = 0; is < 3; is++) { + float* data = reinterpret_cast(buff + nextDynOffs); + LOGP(debug, "splinID={} start offset {} -> {}", is, nextDynOffs, (void*)data); + for (int sector = 0; sector < origCorr.mGeo.getNumberOfSectors(); sector++) { + podMap.mSplineDataOffsets[sector][is] = nextDynOffs; + size_t rowDataOffs = 0; + for (int row = 0; row < NROWS; row++) { + const auto& spline = origCorr.getSpline(sector, row); + const float* dataOr = origCorr.getCorrectionData(sector, row, is); + int nPar = spline.getNumberOfParameters(); + if (is == 1) { + nPar = nPar / 3; + } + if (is == 2) { + nPar = nPar * 2 / 3; + } + LOGP(debug, "Copying {} floats for spline{} of sector:{} row:{} to offset {}", nPar, is, sector, row, nextDynOffs); + size_t nbcopy = nPar * sizeof(float); + if (buffSize < nextDynOffs + nbcopy) { + throw std::runtime_error(fmt::format("attempt to copy {} bytes of data for spline{} of sector{}/row{} to {}, overflowing the buffer of size {}", nbcopy, is, sector, row, nextDynOffs, buffSize)); + } + std::memcpy(data, dataOr, nbcopy); + podMap.getSectorRowInfo(sector, row).dataOffsetBytes[is] = rowDataOffs; + rowDataOffs += nbcopy; + data += nPar; + nextDynOffs += nbcopy; + } + } + } + podMap.mTotalSize = alignOffset(nextDynOffs); + if (buffSize != podMap.mTotalSize) { + throw std::runtime_error(fmt::format("Estimated buffer size {} differs from filled one {}", buffSize, podMap.mTotalSize)); + } + return &getNonConst(buff); +} + +TPCFastTransformPOD* TPCFastTransformPOD::create(char* buff, size_t buffSize, const TPCFastTransform& src) +{ + // instantiate objec to already created buffer of the right size + auto podMap = create(buff, buffSize, src.getCorrection()); + // set data members of TPCFastTransform + podMap->mVdrift = src.getVDrift(); + podMap->mT0 = src.getT0(); + // copy fixed size data --- end + return podMap; +} + +bool TPCFastTransformPOD::test(const TPCFastSpaceChargeCorrection& origCorr, int npoints) const +{ + if (npoints < 1) { + return false; + } + std::vector sector, row; + std::vector y, z; + std::vector> corr0, corr1; + std::vector> corrInv0, corrInv1; + std::vector corrInvX0, corrInvX1; + + sector.reserve(npoints); + row.reserve(npoints); + y.reserve(npoints); + z.reserve(npoints); + corr0.resize(npoints); + corr1.resize(npoints); + corrInv0.resize(npoints); + corrInv1.resize(npoints); + corrInvX0.resize(npoints); + corrInvX1.resize(npoints); + + for (int i = 0; i < npoints; i++) { + sector.push_back(gRandom->Integer(NSECTORS)); + row.push_back(gRandom->Integer(NROWS)); + y.push_back(2 * (gRandom->Rndm() - 0.5) * mGeo.getRowInfo(row.back()).getYmax()); + z.push_back((sector.back() < NSECTORS / 2 ? 1.f : -1.f) * gRandom->Rndm() * 240); + } + long origStart[3], origEnd[3], thisStart[3], thisEnd[3]; + origStart[0] = std::chrono::time_point_cast(std::chrono::system_clock::now()).time_since_epoch().count(); + for (int i = 0; i < npoints; i++) { + corr0.push_back(origCorr.getCorrectionLocal(sector[i], row[i], y[i], z[i])); + } + + origEnd[0] = origStart[1] = std::chrono::time_point_cast(std::chrono::system_clock::now()).time_since_epoch().count(); + for (int i = 0; i < npoints; i++) { + corrInv0.push_back(origCorr.getCorrectionYZatRealYZ(sector[i], row[i], y[i], z[i])); + } + + origEnd[1] = origStart[2] = std::chrono::time_point_cast(std::chrono::system_clock::now()).time_since_epoch().count(); + for (int i = 0; i < npoints; i++) { + corrInvX0.push_back(origCorr.getCorrectionXatRealYZ(sector[i], row[i], y[i], z[i])); + } + // + origEnd[2] = thisStart[0] = std::chrono::time_point_cast(std::chrono::system_clock::now()).time_since_epoch().count(); + for (int i = 0; i < npoints; i++) { + corr1.push_back(this->getCorrectionLocal(sector[i], row[i], y[i], z[i])); + } + thisEnd[0] = thisStart[1] = std::chrono::time_point_cast(std::chrono::system_clock::now()).time_since_epoch().count(); + for (int i = 0; i < npoints; i++) { + corrInv1.push_back(this->getCorrectionYZatRealYZ(sector[i], row[i], y[i], z[i])); + } + + thisEnd[1] = thisStart[2] = std::chrono::time_point_cast(std::chrono::system_clock::now()).time_since_epoch().count(); + for (int i = 0; i < npoints; i++) { + corrInvX1.push_back(this->getCorrectionXatRealYZ(sector[i], row[i], y[i], z[i])); + } + thisEnd[2] = std::chrono::time_point_cast(std::chrono::system_clock::now()).time_since_epoch().count(); + // + size_t ndiff[3] = {}; + for (int i = 0; i < npoints; i++) { + if (corr0[i][0] != corr1[i][0] || corr0[i][1] != corr1[i][1] || corr0[i][2] != corr1[i][2]) { + ndiff[0]++; + } + if (corrInv0[i][0] != corrInv1[i][0] || corrInv0[i][1] != corrInv1[i][1]) { + ndiff[1]++; + } + if (corrInvX0[i] != corrInvX1[i]) { + ndiff[2]++; + } + } + // + LOGP(info, " (ns per call) original this Nmissmatch"); + LOGP(info, "getCorrection {:.3e} {:.3e} {}", double(origEnd[0] - origStart[0]) / npoints * 1000., double(thisEnd[0] - thisStart[0]) / npoints * 1000., ndiff[0]); + LOGP(info, "getCorrectionInvCorrectedX {:.3e} {:.3e} {}", double(origEnd[1] - origStart[1]) / npoints * 1000., double(thisEnd[1] - thisStart[1]) / npoints * 1000., ndiff[1]); + LOGP(info, "getCorrectionInvUV {:.3e} {:.3e} {}", double(origEnd[2] - origStart[2]) / npoints * 1000., double(thisEnd[2] - thisStart[2]) / npoints * 1000., ndiff[2]); + return ndiff[0] == 0 && ndiff[1] == 0 && ndiff[2] == 0; +} + +#endif + +} // namespace gpu +} // namespace o2 diff --git a/GPU/TPCFastTransformation/TPCFastTransformPOD.h b/GPU/TPCFastTransformation/TPCFastTransformPOD.h new file mode 100644 index 0000000000000..0d06a13d0a83c --- /dev/null +++ b/GPU/TPCFastTransformation/TPCFastTransformPOD.h @@ -0,0 +1,920 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file TPCFastTransformPOD.h +/// \brief POD correction map +/// +/// \author ruben.shahoayn@cern.ch + +#ifndef ALICEO2_GPU_TPCFastTransformPOD_H +#define ALICEO2_GPU_TPCFastTransformPOD_H + +#include "TPCFastTransform.h" + +/* +Binary buffer should be cast to TPCFastTransformPOD class using static TPCFastTransformPOD& t = get(buffer); method, +so that its head becomes `this` pointer of the object. + +First we have all the fixed size data members mentioned explicitly. Part of them is duplicating fixed size +data members of TPCFastSpaceChargeCorrection but those starting with mOffs... provide the offset in bytes +(wrt this) for dynamic data which cannot be declared as data member explicitly (since we cannot have any +pointer except `this`) but obtained via getters using stored offsets wrt `this`. +This is followed dynamic part itself. + +dynamic part layout: +1) size_t[ mNumberOfScenarios ] array starting at offset mOffsScenariosOffsets, each element is the offset +of distict spline object (scenario in TPCFastSpaceChargeCorrection) +2) size_t[ mNSplineIDs ] array starting at offset mOffsSplineDataOffsets, each element is the offset of the +beginning of splines data for give splineID + +*/ + +namespace o2 +{ +namespace gpu +{ +class TPCFastTransformPOD +{ + public: + using SliceInfo = TPCFastSpaceChargeCorrection::SliceInfo; // obsolete + using GridInfo = TPCFastSpaceChargeCorrection::GridInfo; + using SectorRowInfo = TPCFastSpaceChargeCorrection::SectorRowInfo; + + using SplineTypeXYZ = TPCFastSpaceChargeCorrection::SplineTypeXYZ; + using SplineTypeInvX = TPCFastSpaceChargeCorrection::SplineTypeInvX; + using SplineTypeInvYZ = TPCFastSpaceChargeCorrection::SplineTypeInvYZ; + using SplineType = TPCFastSpaceChargeCorrection::SplineType; + + /// convert prefilled buffer to TPCFastTransformPOD + GPUd() static const TPCFastTransformPOD& get(const char* head) { return *reinterpret_cast(head); } + + /// _______________ high level methods a la TPCFastTransform _______________________ + /// + // Methods taking extra reference transform are legacy compound transforms used to scale corrections. + GPUd() void Transform(int32_t sector, int32_t row, float pad, float time, float& x, float& y, float& z, float vertexTime = 0, const TPCFastTransformPOD* ref = nullptr, const TPCFastTransformPOD* ref2 = nullptr, float scale = 0.f, float scale2 = 0.f, int32_t scaleMode = 0) const; + GPUd() void TransformXYZ(int32_t sector, int32_t row, float& x, float& y, float& z, const TPCFastTransformPOD* ref = nullptr, const TPCFastTransformPOD* ref2 = nullptr, float scale = 0.f, float scale2 = 0.f, int32_t scaleMode = 0) const; + + GPUd() void Transform_new(int32_t sector, int32_t row, float pad, float time, float& x, float& y, float& z, float vertexTime = 0) const; + GPUd() void TransformXYZ_new(int32_t sector, int32_t row, float& x, float& y, float& z) const; + + /// Transformation in the time frame + GPUd() void TransformInTimeFrame(int32_t sector, int32_t row, float pad, float time, float& x, float& y, float& z, float maxTimeBin) const; + GPUd() void TransformInTimeFrame(int32_t sector, float time, float& z, float maxTimeBin) const; + + /// Inverse transformation + GPUd() void InverseTransformInTimeFrame(int32_t sector, int32_t row, float /*x*/, float y, float z, float& pad, float& time, float maxTimeBin) const; + GPUd() float InverseTransformInTimeFrame(int32_t sector, float z, float maxTimeBin) const; + + /// Inverse transformation: Transformed Y and Z -> transformed X + GPUd() void InverseTransformYZtoX(int32_t sector, int32_t row, float y, float z, float& x, const TPCFastTransformPOD* ref = nullptr, const TPCFastTransformPOD* ref2 = nullptr, float scale = 0.f, float scale2 = 0.f, int32_t scaleMode = 0) const; + GPUd() void InverseTransformYZtoX_new(int32_t sector, int32_t row, float y, float z, float& x) const; + + /// Inverse transformation: Transformed Y and Z -> Y and Z, transformed w/o space charge correction + GPUd() void InverseTransformYZtoNominalYZ(int32_t sector, int32_t row, float y, float z, float& ny, float& nz, const TPCFastTransformPOD* ref = nullptr, const TPCFastTransformPOD* ref2 = nullptr, float scale = 0.f, float scale2 = 0.f, int32_t scaleMode = 0) const; + GPUd() void InverseTransformYZtoNominalYZ_new(int32_t sector, int32_t row, float y, float z, float& ny, float& nz) const; + + /// Inverse transformation: Transformed X, Y and Z -> X, Y and Z, transformed w/o space charge correction + GPUd() void InverseTransformXYZtoNominalXYZ(int32_t sector, int32_t row, float x, float y, float z, float& nx, float& ny, float& nz, const TPCFastTransformPOD* ref = nullptr, const TPCFastTransformPOD* ref2 = nullptr, float scale = 0.f, float scale2 = 0.f, int32_t scaleMode = 0) const; + GPUd() void InverseTransformXYZtoNominalXYZ_new(int32_t sector, int32_t row, float x, float y, float z, float& nx, float& ny, float& nz) const; + + /// Ideal transformation with Vdrift only - without calibration + GPUd() void TransformIdeal(int32_t sector, int32_t row, float pad, float time, float& x, float& y, float& z, float vertexTime) const; + GPUd() void TransformIdealZ(int32_t sector, float time, float& z, float vertexTime) const; + + GPUd() void convPadTimeToLocal(int32_t sector, int32_t row, float pad, float time, float& y, float& z, float vertexTime) const; + GPUd() void convPadTimeToLocalInTimeFrame(int32_t sector, int32_t row, float pad, float time, float& y, float& z, float maxTimeBin) const; + + GPUd() void convLocalToPadTime(int32_t sector, int32_t row, float y, float z, float& pad, float& time, float vertexTime) const; + GPUd() void convLocalToPadTimeInTimeFrame(int32_t sector, int32_t row, float y, float z, float& pad, float& time, float maxTimeBin) const; + + GPUd() float convTimeToZinTimeFrame(int32_t sector, float time, float maxTimeBin) const; + GPUd() float convZtoTimeInTimeFrame(int32_t sector, float z, float maxTimeBin) const; + GPUd() float convDeltaTimeToDeltaZinTimeFrame(int32_t sector, float deltaTime) const; + GPUd() float convDeltaZtoDeltaTimeInTimeFrame(int32_t sector, float deltaZ) const; + GPUd() float convDeltaZtoDeltaTimeInTimeFrameAbs(float deltaZ) const; + GPUd() float convZOffsetToVertexTime(int32_t sector, float zOffset, float maxTimeBin) const; + GPUd() float convVertexTimeToZOffset(int32_t sector, float vertexTime, float maxTimeBin) const; + + /// _______________ methods a la TPCFastSpaceChargeCorrection: cluster correction _______________________ + void setApplyCorrectionOn() { mApplyCorrection = 1; } + void setApplyCorrectionOff() { mApplyCorrection = 0; } + bool isCorrectionApplied() { return mApplyCorrection; } + + /// TPC geometry information + GPUd() const TPCFastTransformGeo& getGeometry() const { return mGeo; } + + /// Gives TPC sector & row info + GPUd() const SectorRowInfo& getSectorRowInfo(int32_t sector, int32_t row) const { return mSectorRowInfos[NROWS * sector + row]; } + + /// Gives TPC sector & row info + GPUd() SectorRowInfo& getSectorRowInfo(int32_t sector, int32_t row) { return mSectorRowInfos[NROWS * sector + row]; } + + /// Gives its own size including dynamic part + GPUd() size_t size() const { return mTotalSize; } + + /// Gives the time stamp of the current calibaration parameters + GPUd() long int getTimeStamp() const { return mTimeStamp; } + + /// Return mVDrift in cm / time bin + GPUd() float getVDrift() const { return mVdrift; } + + /// Return T0 in time bin units + GPUd() float getT0() const { return mT0; } + + /// Return IDC estimator + GPUd() float getIDC() const { return mIDC; } + + /// Return Lumi estimator + GPUd() float getLumi() const { return mLumi; } + + /// maximal possible drift time of the active area + GPUd() float getMaxDriftTime(int32_t sector, int32_t row, float pad) const; + + /// maximal possible drift time of the active area + GPUd() float getMaxDriftTime(int32_t sector, int32_t row) const; + + /// maximal possible drift time of the active area + GPUd() float getMaxDriftTime(int32_t sector) const; + + /// Sets the time stamp of the current calibaration + GPUd() void setTimeStamp(long int v) { mTimeStamp = v; } + + /// Sets current vdrift + GPUd() void setVDrift(float v) { mVdrift = v; } + + /// Sets current T0 + GPUd() void setT0(float v) { mT0 = v; } + + /// Sets IDC estimator + GPUd() void setIDC(float v) { mIDC = v; } + + /// Sets CTP Lumi estimator + GPUd() void setLumi(float v) { mLumi = v; } + + /// Gives a reference to a spline + GPUd() const SplineType& getSpline(int32_t sector, int32_t row) const { return *reinterpret_cast(getThis() + getScenarioOffset(getSectorRowInfo(sector, row).splineScenarioID)); } + + /// Gives pointer to spline data + GPUd() const float* getCorrectionData(int32_t sector, int32_t row, int32_t iSpline = 0) const { return reinterpret_cast(getThis() + mSplineDataOffsets[sector][iSpline] + getSectorRowInfo(sector, row).dataOffsetBytes[iSpline]); } + + /// Gives const pointer to a spline for the inverse X correction + GPUd() const SplineTypeInvX& getSplineInvX(int32_t sector, int32_t row) const { return reinterpret_cast(getSpline(sector, row)); } + + /// Gives pointer to spline data for the inverse X correction + GPUd() const float* getCorrectionDataInvX(int32_t sector, int32_t row) const { return getCorrectionData(sector, row, 1); } + + /// Gives const pointer to a spline for the inverse YZ correction + GPUd() const SplineTypeInvYZ& getSplineInvYZ(int32_t sector, int32_t row) const { return reinterpret_cast(getSpline(sector, row)); } + + /// Gives pointer to spline data for the inverse YZ correction + GPUd() const float* getCorrectionDataInvYZ(int32_t sector, int32_t row) const { return getCorrectionData(sector, row, 2); } + + /// _______________ The main method: cluster correction _______________________ + GPUdi() std::array getCorrectionLocal(int32_t sector, int32_t row, float y, float z) const; + + /// inverse correction: Real Y and Z -> Real X + GPUd() float getCorrectionXatRealYZ(int32_t sector, int32_t row, float realY, float realZ) const; + + /// inverse correction: Real Y and Z -> measred Y and Z + GPUd() std::array getCorrectionYZatRealYZ(int32_t sector, int32_t row, float realY, float realZ) const; + + /// transformation in the sector local frame + GPUd() void TransformLocal(int32_t sector, int32_t row, float& x, float& y, float& z, const TPCFastTransformPOD* ref, const TPCFastTransformPOD* ref2, float scale, float scale2, int32_t scaleMode) const; + GPUd() void TransformLocal(int32_t sector, int32_t row, float& x, float& y, float& z) const; + + /// _______________ Utilities _______________________________________________ + + /// convert local y, z to internal grid coordinates u,v + /// return values: u, v, scaling factor + GPUd() std::array convLocalToGrid(int32_t sector, int32_t row, float y, float z) const; + + /// convert internal grid coordinates u,v to local y, z + /// return values: y, z, scaling factor + GPUd() std::array convGridToLocal(int32_t sector, int32_t row, float u, float v) const; + + /// convert real Y, Z to the internal grid coordinates + /// return values: u, v, scaling factor + GPUd() std::array convRealLocalToGrid(int32_t sector, int32_t row, float y, float z) const; + + /// convert internal grid coordinates to the real Y, Z + /// return values: y, z + GPUd() std::array convGridToRealLocal(int32_t sector, int32_t row, float u, float v) const; + + GPUd() bool isLocalInsideGrid(int32_t sector, int32_t row, float y, float z) const; + GPUd() bool isRealLocalInsideGrid(int32_t sector, int32_t row, float y, float z) const; + +#if !defined(GPUCA_GPUCODE) + /// Create POD transform from old flat-buffer one. Provided vector will serve as a buffer + template + static TPCFastTransformPOD* create(V& destVector, const TPCFastTransform& src); + + /// create filling only part corresponding to TPCFastSpaceChargeCorrection. Data members coming from TPCFastTransform (e.g. VDrift, T0..) are not set + template + static TPCFastTransformPOD* create(V& destVector, const TPCFastSpaceChargeCorrection& src); + + bool test(const TPCFastTransform& src, int32_t npoints = 100000) const { return test(src.getCorrection(), npoints); } + bool test(const TPCFastSpaceChargeCorrection& origCorr, int32_t npoints = 100000) const; +#endif + + /// Print method + void print() const; + + GPUd() float convDriftLengthToTime(float driftLength, float vertexTime) const; + + static constexpr int NROWS = 152; + static constexpr int NSECTORS = TPCFastTransformGeo::getNumberOfSectors(); + static constexpr int NSplineIDs = 3; ///< number of spline data sets for each sector/row + + private: +#if !defined(GPUCA_GPUCODE) + static constexpr size_t AlignmentBytes = 8; + static size_t alignOffset(size_t offs) + { + auto res = offs % AlignmentBytes; + return res ? offs + (AlignmentBytes - res) : offs; + } + static size_t estimateSize(const TPCFastTransform& src) { return estimateSize(src.getCorrection()); } + static size_t estimateSize(const TPCFastSpaceChargeCorrection& origCorr); + static TPCFastTransformPOD* create(char* buff, size_t buffSize, const TPCFastTransform& src); + static TPCFastTransformPOD* create(char* buff, size_t buffSize, const TPCFastSpaceChargeCorrection& src); + ///< get address to which the offset in bytes must be added to arrive to particular dynamic part + GPUd() const char* getThis() const { return reinterpret_cast(this); } + GPUd() static TPCFastTransformPOD& getNonConst(char* head) { return *reinterpret_cast(head); } +#endif + + ///< return offset of the spline object start (equivalent of mScenarioPtr in the TPCFastSpaceChargeCorrection) + GPUd() const size_t getScenarioOffset(int s) const { return (reinterpret_cast(getThis() + mOffsScenariosOffsets))[s]; } + + bool mApplyCorrection{}; ///< flag to apply corrections + int mNumberOfScenarios{}; ///< Number of approximation spline scenarios + size_t mTotalSize{}; ///< total size of the buffer + size_t mOffsScenariosOffsets{}; ///< start of the array of mNumberOfScenarios offsets for each type of spline + size_t mSplineDataOffsets[TPCFastTransformGeo::getNumberOfSectors()][NSplineIDs]; ///< start of data for each sector and iSpline data + long int mTimeStamp{}; ///< time stamp of the current calibration + float mT0; ///< T0 in [time bin] + float mVdrift; ///< VDrift in [cm/time bin] + float mLumi; ///< luminosity estimator (for info only) + float mIDC; ///< IDC estimator (for info only) + + TPCFastTransformGeo mGeo; ///< TPC geometry information + SectorRowInfo mSectorRowInfos[NROWS * TPCFastTransformGeo::getNumberOfSectors()]; + + ClassDefNV(TPCFastTransformPOD, 0); +}; + +GPUdi() std::array TPCFastTransformPOD::getCorrectionLocal(int32_t sector, int32_t row, float y, float z) const +{ + const auto& info = getSectorRowInfo(sector, row); + const SplineType& spline = getSpline(sector, row); + const float* splineData = getCorrectionData(sector, row); + + auto val = convLocalToGrid(sector, row, y, z); + + float dxyz[3]; + spline.interpolateAtU(splineData, val[0], val[1], dxyz); + + if (CAMath::Abs(dxyz[0]) > 100.f || CAMath::Abs(dxyz[1]) > 100.f || CAMath::Abs(dxyz[2]) > 100.f) { + val[2] = 0.f; // TODO: DR: Protect from FPEs, fix upstream and remove once guaranteed that it is fixed + } + + float dx = val[2] * GPUCommonMath::Clamp(dxyz[0], info.minCorr[0], info.maxCorr[0]); + float dy = val[2] * GPUCommonMath::Clamp(dxyz[1], info.minCorr[1], info.maxCorr[1]); + float dz = val[2] * GPUCommonMath::Clamp(dxyz[2], info.minCorr[2], info.maxCorr[2]); + + return {dx, dy, dz}; +} + +GPUdi() float TPCFastTransformPOD::getCorrectionXatRealYZ(int32_t sector, int32_t row, float realY, float realZ) const +{ + const auto& info = getSectorRowInfo(sector, row); + auto val = convRealLocalToGrid(sector, row, realY, realZ); + float dx = 0; + getSplineInvX(sector, row).interpolateAtU(getCorrectionDataInvX(sector, row), val[0], val[1], &dx); + if (CAMath::Abs(dx) > 100.f) { + val[2] = 0.f; // TODO: DR: Protect from FPEs, fix upstream and remove once guaranteed that it is fixed + } + dx = val[2] * GPUCommonMath::Clamp(dx, info.minCorr[0], info.maxCorr[0]); + return dx; +} + +GPUdi() std::array TPCFastTransformPOD::getCorrectionYZatRealYZ(int32_t sector, int32_t row, float realY, float realZ) const +{ + auto val = convRealLocalToGrid(sector, row, realY, realZ); + const auto& info = getSectorRowInfo(sector, row); + float dyz[2]; + getSplineInvYZ(sector, row).interpolateAtU(getCorrectionDataInvYZ(sector, row), val[0], val[1], dyz); + if (CAMath::Abs(dyz[0]) > 100.f || CAMath::Abs(dyz[1]) > 100.f) { + val[2] = 0.f; // TODO: DR: Protect from FPEs, fix upstream and remove once guaranteed that it is fixed + } + dyz[0] = val[2] * GPUCommonMath::Clamp(dyz[0], info.minCorr[1], info.maxCorr[1]); + dyz[1] = val[2] * GPUCommonMath::Clamp(dyz[1], info.minCorr[2], info.maxCorr[2]); + return {dyz[0], dyz[1]}; +} + +GPUdi() std::array TPCFastTransformPOD::convLocalToGrid(int32_t sector, int32_t row, float y, float z) const +{ + /// convert local y, z to internal grid coordinates u,v + /// return values: u, v, scaling factor + const SplineType& spline = getSpline(sector, row); + auto val = getSectorRowInfo(sector, row).gridMeasured.convLocalToGridUntruncated(y, z); + // shrink to the grid + val[0] = GPUCommonMath::Clamp(val[0], 0.f, (float)spline.getGridX1().getUmax()); + val[1] = GPUCommonMath::Clamp(val[1], 0.f, (float)spline.getGridX2().getUmax()); + return val; +} + +GPUdi() std::array TPCFastTransformPOD::convGridToLocal(int32_t sector, int32_t row, float gridU, float gridV) const +{ + /// convert internal grid coordinates u,v to local y, z + return getSectorRowInfo(sector, row).gridMeasured.convGridToLocal(gridU, gridV); +} + +GPUdi() std::array TPCFastTransformPOD::convRealLocalToGrid(int32_t sector, int32_t row, float y, float z) const +{ + /// convert real y, z to the internal grid coordinates + scale + const SplineType& spline = getSpline(sector, row); + auto val = getSectorRowInfo(sector, row).gridReal.convLocalToGridUntruncated(y, z); + // shrink to the grid + val[0] = GPUCommonMath::Clamp(val[0], 0.f, (float)spline.getGridX1().getUmax()); + val[1] = GPUCommonMath::Clamp(val[1], 0.f, (float)spline.getGridX2().getUmax()); + return val; +} + +GPUdi() std::array TPCFastTransformPOD::convGridToRealLocal(int32_t sector, int32_t row, float gridU, float gridV) const +{ + /// convert internal grid coordinates u,v to the real y, z + return getSectorRowInfo(sector, row).gridReal.convGridToLocal(gridU, gridV); +} + +GPUdi() bool TPCFastTransformPOD::isLocalInsideGrid(int32_t sector, int32_t row, float y, float z) const +{ + /// check if local y, z are inside the grid + auto val = getSectorRowInfo(sector, row).gridMeasured.convLocalToGridUntruncated(y, z); + const auto& spline = getSpline(sector, row); + // shrink to the grid + if (val[0] < 0.f || val[0] > (float)spline.getGridX1().getUmax() || // + val[1] < 0.f || val[1] > (float)spline.getGridX2().getUmax()) { + return false; + } + return true; +} + +GPUdi() bool TPCFastTransformPOD::isRealLocalInsideGrid(int32_t sector, int32_t row, float y, float z) const +{ + /// check if local y, z are inside the grid + auto val = getSectorRowInfo(sector, row).gridReal.convLocalToGridUntruncated(y, z); + const auto& spline = getSpline(sector, row); + // shrink to the grid + if (val[0] < 0.f || val[0] > (float)spline.getGridX1().getUmax() || // + val[1] < 0.f || val[1] > (float)spline.getGridX2().getUmax()) { + return false; + } + return true; +} + +#if !defined(GPUCA_GPUCODE) +/// Create POD transform from old flat-buffer one. Provided vector will serve as a buffer +template +TPCFastTransformPOD* TPCFastTransformPOD::create(V& destVector, const TPCFastTransform& src) +{ + const auto& origCorr = src.getCorrection(); + size_t estSize = estimateSize(src); + destVector.resize(estSize); // allocate exact size + LOGP(debug, "OrigCorrSize:{} SelfSize: {} Estimated POS size: {}", src.getCorrection().getFlatBufferSize(), sizeof(TPCFastTransformPOD), estSize); + char* base = destVector.data(); + auto res = create(destVector.data(), destVector.size(), src); + res->setTimeStamp(src.getTimeStamp()); + res->setVDrift(src.getVDrift()); + res->setT0(src.getT0()); + res->setLumi(src.getLumi()); + res->setIDC(src.getIDC()); + return res; +} + +template +TPCFastTransformPOD* TPCFastTransformPOD::create(V& destVector, const TPCFastSpaceChargeCorrection& origCorr) +{ + // create filling only part corresponding to TPCFastSpaceChargeCorrection. Data members coming from TPCFastTransform (e.g. VDrift, T0..) are not set + size_t estSize = estimateSize(origCorr); + destVector.resize(estSize); // allocate exact size + LOGP(debug, "OrigCorrSize:{} SelfSize: {} Estimated POS size: {}", origCorr.getFlatBufferSize(), sizeof(TPCFastTransformPOD), estSize); + char* base = destVector.data(); + return create(destVector.data(), destVector.size(), origCorr); +} +#endif + +GPUdi() void TPCFastTransformPOD::TransformLocal(int32_t sector, int32_t row, float& x, float& y, float& z, const TPCFastTransformPOD* ref, const TPCFastTransformPOD* ref2, float scale, float scale2, int32_t scaleMode) const +{ + GPUCA_RTC_SPECIAL_CODE(ref2 = nullptr; scale2 = 0.f;); + + if (!mApplyCorrection) { + return; + } + + float dx = 0.f, dy = 0.f, dz = 0.f; + + if ((scale >= 0.f) || (scaleMode == 1) || (scaleMode == 2)) { + const auto corrLocal = getCorrectionLocal(sector, row, y, z); + dx = corrLocal[0]; + dy = corrLocal[1]; + dz = corrLocal[2]; + if (ref) { + if ((scale > 0.f) && (scaleMode == 0)) { // scaling was requested + auto val = ref->getCorrectionLocal(sector, row, y, z); + dx = (dx - val[0]) * scale + val[0]; + dy = (dy - val[1]) * scale + val[1]; + dz = (dz - val[2]) * scale + val[2]; + } else if ((scale != 0.f) && ((scaleMode == 1) || (scaleMode == 2))) { + auto val = ref->getCorrectionLocal(sector, row, y, z); + dx = val[0] * scale + dx; + dy = val[1] * scale + dy; + dz = val[2] * scale + dz; + } + } + if (ref2 && (scale2 != 0)) { + auto val = ref2->getCorrectionLocal(sector, row, y, z); + dx = val[0] * scale2 + dx; + dy = val[1] * scale2 + dy; + dz = val[2] * scale2 + dz; + } + } + + GPUCA_DEBUG_STREAMER_CHECK(if (o2::utils::DebugStreamer::checkStream(o2::utils::StreamFlags::streamFastTransform)) { + float lx = x, ly = y, lz = z; + + float gx, gy, gz; + getGeometry().convLocalToGlobal(sector, lx, ly, lz, gx, gy, gz); + + float lxT = lx + dx; + float lyT = ly + dy; + float lzT = lz + dz; + + float invYZtoXScaled; + InverseTransformYZtoX(sector, row, lyT, lzT, invYZtoXScaled, ref, ref2, scale, scale2, scaleMode); + + float invYZtoX; + InverseTransformYZtoX(sector, row, lyT, lzT, invYZtoX); + + float YZtoNominalY; + float YZtoNominalZ; + InverseTransformYZtoNominalYZ(sector, row, lyT, lzT, YZtoNominalY, YZtoNominalZ); + + float YZtoNominalYScaled; + float YZtoNominalZScaled; + InverseTransformYZtoNominalYZ(sector, row, lyT, lzT, YZtoNominalYScaled, YZtoNominalZScaled, ref, ref2, scale, scale2, scaleMode); + + float dxRef = 0.f, dyRef = 0.f, dzRef = 0.f; + if (ref) { + const auto corr = ref->getCorrectionLocal(sector, row, y, z); + dxRef = corr[0]; + dyRef = corr[1]; + dzRef = corr[2]; + } + + float dxRef2 = 0.f, dyRef2 = 0.f, dzRef2 = 0.f; + if (ref2) { + const auto corr = ref2->getCorrectionLocal(sector, row, y, z); + dxRef2 = corr[0]; + dyRef2 = corr[1]; + dzRef2 = corr[2]; + } + + auto [dxOrig, dyOrig, dzOrig] = getCorrectionLocal(sector, row, y, z); + + o2::utils::DebugStreamer::instance()->getStreamer("debug_fasttransform", "UPDATE") << o2::utils::DebugStreamer::instance()->getUniqueTreeName("tree_Transform").data() + // corrections in x, u, v + << "dxOrig=" << dxOrig + << "dyOrig=" << dyOrig + << "dzOrig=" << dzOrig + << "dxRef=" << dxRef + << "dyRef=" << dyRef + << "dzRef=" << dzRef + << "dxRef2=" << dxRef2 + << "dyRef2=" << dyRef2 + << "dzRef2=" << dzRef2 + << "dx=" << dx + << "dy=" << dy + << "dz=" << dz + << "row=" << row + << "sector=" << sector + << "scale=" << scale + << "scale2=" << scale2 + // original local coordinates + << "ly=" << ly + << "lz=" << lz + << "lx=" << lx + // corrected local coordinated + << "lxT=" << lxT + << "lyT=" << lyT + << "lzT=" << lzT + // global uncorrected coordinates + << "gx=" << gx + << "gy=" << gy + << "gz=" << gz + // some transformations which are applied + << "invYZtoX=" << invYZtoX + << "invYZtoXScaled=" << invYZtoXScaled + << "YZtoNominalY=" << YZtoNominalY + << "YZtoNominalYScaled=" << YZtoNominalYScaled + << "YZtoNominalZ=" << YZtoNominalZ + << "YZtoNominalZScaled=" << YZtoNominalZScaled + << "scaleMode=" << scaleMode + << "\n"; + }) + + x += dx; + y += dy; + z += dz; +} + +GPUdi() void TPCFastTransformPOD::TransformLocal(int32_t sector, int32_t row, float& x, float& y, float& z) const +{ + if (!mApplyCorrection) { + return; + } + const auto corrLocal = getCorrectionLocal(sector, row, y, z); + + GPUCA_DEBUG_STREAMER_CHECK(if (o2::utils::DebugStreamer::checkStream(o2::utils::StreamFlags::streamFastTransform)) { + float lx = x, ly = y, lz = z; + float gx, gy, gz; + getGeometry().convLocalToGlobal(sector, lx, ly, lz, gx, gy, gz); + dx = corrLocal[0]; + dy = corrLocal[1]; + dz = corrLocal[2]; + float lxT = lx + dx; + float lyT = ly + dy; + float lzT = lz + dz; + float invYZtoX; + InverseTransformYZtoX_new(sector, row, lyT, lzT, invYZtoX); + + float YZtoNominalY; + float YZtoNominalZ; + InverseTransformYZtoNominalYZ_new(sector, row, lyT, lzT, YZtoNominalY, YZtoNominalZ); + + o2::utils::DebugStreamer::instance()->getStreamer("debug_fasttransform", "UPDATE") << o2::utils::DebugStreamer::instance()->getUniqueTreeName("tree_Transform").data() + // corrections in x, u, v + << "dx=" << dx + << "dy=" << dy + << "dz=" << dz + << "row=" << row + << "sector=" << sector + // original local coordinates + << "ly=" << ly + << "lz=" << lz + << "lx=" << lx + // corrected local coordinated + << "lxT=" << lxT + << "lyT=" << lyT + << "lzT=" << lzT + // global uncorrected coordinates + << "gx=" << gx + << "gy=" << gy + << "gz=" << gz + // some transformations which are applied + << "invYZtoX=" << invYZtoX + << "YZtoNominalY=" << YZtoNominalY + << "YZtoNominalZ=" << YZtoNominalZ + << "\n"; + }) + + x += corrLocal[0]; + y += corrLocal[1]; + z += corrLocal[2]; +} + +GPUdi() void TPCFastTransformPOD::Transform(int32_t sector, int32_t row, float pad, float time, float& x, float& y, float& z, float vertexTime, const TPCFastTransformPOD* ref, const TPCFastTransformPOD* ref2, float scale, float scale2, int32_t scaleMode) const +{ + /// _______________ The main method: cluster transformation _______________________ + /// + /// Transforms raw TPC coordinates to local XYZ withing a sector + /// taking calibration into account. + /// + + const TPCFastTransformGeo::RowInfo& rowInfo = getGeometry().getRowInfo(row); + + x = rowInfo.x; + convPadTimeToLocal(sector, row, pad, time, y, z, vertexTime); + TransformLocal(sector, row, x, y, z, ref, ref2, scale, scale2, scaleMode); +} + +GPUdi() void TPCFastTransformPOD::Transform_new(int32_t sector, int32_t row, float pad, float time, float& x, float& y, float& z, float vertexTime) const +{ + /// _______________ The main method: cluster transformation _______________________ + /// + /// Transforms raw TPC coordinates to local XYZ withing a sector + /// taking calibration into account. + /// + + const TPCFastTransformGeo::RowInfo& rowInfo = getGeometry().getRowInfo(row); + + x = rowInfo.x; + convPadTimeToLocal(sector, row, pad, time, y, z, vertexTime); + TransformLocal(sector, row, x, y, z); +} + +GPUdi() void TPCFastTransformPOD::TransformXYZ(int32_t sector, int32_t row, float& x, float& y, float& z, const TPCFastTransformPOD* ref, const TPCFastTransformPOD* ref2, float scale, float scale2, int32_t scaleMode) const +{ + + TransformLocal(sector, row, x, y, z, ref, ref2, scale, scale2, scaleMode); +} + +GPUdi() void TPCFastTransformPOD::TransformXYZ_new(int32_t sector, int32_t row, float& x, float& y, float& z) const +{ + + TransformLocal(sector, row, x, y, z); +} + +GPUdi() void TPCFastTransformPOD::TransformInTimeFrame(int32_t sector, float time, float& z, float maxTimeBin) const +{ + float l = (time - mT0 - maxTimeBin) * mVdrift; // drift length cm + z = getGeometry().convDriftLengthToZ1(sector, l); +} + +GPUdi() void TPCFastTransformPOD::TransformInTimeFrame(int32_t sector, int32_t row, float pad, float time, float& x, float& y, float& z, float maxTimeBin) const +{ + /// _______________ Special cluster transformation for a time frame _______________________ + /// + /// Same as Transform(), but clusters are shifted in z such, that Z(maxTimeBin)==0 + /// Corrections and Time-Of-Flight correction are not alpplied. + /// + + const TPCFastTransformGeo::RowInfo& rowInfo = getGeometry().getRowInfo(row); + x = rowInfo.x; + convPadTimeToLocalInTimeFrame(sector, row, pad, time, y, z, maxTimeBin); +} + +GPUdi() void TPCFastTransformPOD::InverseTransformInTimeFrame(int32_t sector, int32_t row, float /*x*/, float y, float z, float& pad, float& time, float maxTimeBin) const +{ + /// Inverse transformation to TransformInTimeFrame + convLocalToPadTimeInTimeFrame(sector, row, y, z, pad, time, maxTimeBin); +} + +GPUdi() float TPCFastTransformPOD::InverseTransformInTimeFrame(int32_t sector, float z, float maxTimeBin) const +{ + float pad, time; + InverseTransformInTimeFrame(sector, 0, 0, 0, z, pad, time, maxTimeBin); + return time; +} + +GPUdi() void TPCFastTransformPOD::TransformIdealZ(int32_t sector, float time, float& z, float vertexTime) const +{ + /// _______________ The main method: cluster transformation _______________________ + /// + /// Transforms time TPC coordinates to local Z withing a sector + /// Ideal transformation: only Vdrift from DCS. + /// No space charge corrections, no time of flight correction + /// + + float l = (time - mT0 - vertexTime) * mVdrift; // drift length cm + z = getGeometry().convDriftLengthToZ1(sector, l); +} + +GPUdi() void TPCFastTransformPOD::TransformIdeal(int32_t sector, int32_t row, float pad, float time, float& x, float& y, float& z, float vertexTime) const +{ + /// _______________ The main method: cluster transformation _______________________ + /// + /// Transforms raw TPC coordinates to local XYZ withing a sector + /// Ideal transformation: only Vdrift from DCS. + /// No space charge corrections, no time of flight correction + /// + + x = getGeometry().getRowInfo(row).x; + float driftLength = (time - mT0 - vertexTime) * mVdrift; // drift length cm + const auto localval = getGeometry().convPadDriftLengthToLocal(sector, row, pad, driftLength); + y = localval[0]; + z = localval[1]; +} + +GPUdi() float TPCFastTransformPOD::convTimeToZinTimeFrame(int32_t sector, float time, float maxTimeBin) const +{ + /// _______________ Special cluster transformation for a time frame _______________________ + /// + /// Same as Transform(), but clusters are shifted in z such, that Z(maxTimeBin)==0 + /// Corrections and Time-Of-Flight correction are not alpplied. + /// Only Z coordinate. + /// + + float v = (time - mT0 - maxTimeBin) * mVdrift; // drift length cm + float z = (sector < getGeometry().getNumberOfSectorsA()) ? -v : v; + return z; +} + +GPUdi() float TPCFastTransformPOD::convZtoTimeInTimeFrame(int32_t sector, float z, float maxTimeBin) const +{ + /// Inverse transformation of convTimeToZinTimeFrame() + float v = (sector < getGeometry().getNumberOfSectorsA()) ? -z : z; + return mT0 + maxTimeBin + v / mVdrift; +} + +GPUdi() float TPCFastTransformPOD::convDeltaTimeToDeltaZinTimeFrame(int32_t sector, float deltaTime) const +{ + float deltaZ = deltaTime * mVdrift; + return sector < getGeometry().getNumberOfSectorsA() ? -deltaZ : deltaZ; +} + +GPUdi() float TPCFastTransformPOD::convDeltaZtoDeltaTimeInTimeFrameAbs(float deltaZ) const +{ + return deltaZ / mVdrift; +} + +GPUdi() float TPCFastTransformPOD::convDeltaZtoDeltaTimeInTimeFrame(int32_t sector, float deltaZ) const +{ + float deltaT = deltaZ / mVdrift; + return sector < getGeometry().getNumberOfSectorsA() ? -deltaT : deltaT; +} + +GPUdi() float TPCFastTransformPOD::getMaxDriftTime(int32_t sector, int32_t row, float pad) const +{ + /// maximal possible drift time of the active area + return convDriftLengthToTime(getGeometry().getTPCzLength(), 0.f); +} + +GPUdi() float TPCFastTransformPOD::getMaxDriftTime(int32_t sector, int32_t row) const +{ + /// maximal possible drift time of the active area + return convDriftLengthToTime(getGeometry().getTPCzLength(), 0.f); +} + +GPUdi() float TPCFastTransformPOD::getMaxDriftTime(int32_t sector) const +{ + /// maximal possible drift time of the active area + return convDriftLengthToTime(getGeometry().getTPCzLength(), 0.f); +} + +GPUdi() void TPCFastTransformPOD::InverseTransformYZtoX(int32_t sector, int32_t row, float realY, float realZ, float& realX, const TPCFastTransformPOD* ref, const TPCFastTransformPOD* ref2, float scale, float scale2, int32_t scaleMode) const +{ + GPUCA_RTC_SPECIAL_CODE(ref2 = nullptr; scale2 = 0.f;); + /// Transformation y,z -> x + + float dx = 0.f; + + if ((scale >= 0.f) || (scaleMode == 1) || (scaleMode == 2)) { + dx = getCorrectionXatRealYZ(sector, row, realY, realZ); + if (ref) { // scaling was requested + if (scaleMode == 0 && scale > 0.f) { + float dxref = ref->getCorrectionXatRealYZ(sector, row, realY, realZ); + dx = (dx - dxref) * scale + dxref; + } else if ((scale != 0) && ((scaleMode == 1) || (scaleMode == 2))) { + float dxref = ref->getCorrectionXatRealYZ(sector, row, realY, realZ); + dx = dxref * scale + dx; + } + } + if (ref2 && (scale2 != 0)) { + float dxref = ref2->getCorrectionXatRealYZ(sector, row, realY, realZ); + dx = dxref * scale2 + dx; + } + } + + realX = getGeometry().getRowInfo(row).x + dx; + + GPUCA_DEBUG_STREAMER_CHECK(if (o2::utils::DebugStreamer::checkStream(o2::utils::StreamFlags::streamFastTransform)) { + o2::utils::DebugStreamer::instance()->getStreamer("debug_fasttransform", "UPDATE") << o2::utils::DebugStreamer::instance()->getUniqueTreeName("tree_InverseTransformYZtoX").data() + << "sector=" << sector + << "row=" << row + << "scale=" << scale + << "y=" << realY + << "z=" << realZ + << "x=" << realX + << "\n"; + }) +} + +GPUdi() void TPCFastTransformPOD::InverseTransformYZtoX_new(int32_t sector, int32_t row, float realY, float realZ, float& realX) const +{ + /// Transformation y,z -> x + + float dx = 0.f; + + dx = getCorrectionXatRealYZ(sector, row, realY, realZ); + realX = getGeometry().getRowInfo(row).x + dx; + + GPUCA_DEBUG_STREAMER_CHECK(if (o2::utils::DebugStreamer::checkStream(o2::utils::StreamFlags::streamFastTransform)) { + o2::utils::DebugStreamer::instance()->getStreamer("debug_fasttransform", "UPDATE") << o2::utils::DebugStreamer::instance()->getUniqueTreeName("tree_InverseTransformYZtoX").data() + << "sector=" << sector + << "row=" << row + << "y=" << realY + << "z=" << realZ + << "x=" << realX + << "\n"; + }) +} + +GPUdi() void TPCFastTransformPOD::InverseTransformYZtoNominalYZ(int32_t sector, int32_t row, float realY, float realZ, float& measuredY, float& measuredZ, const TPCFastTransformPOD* ref, const TPCFastTransformPOD* ref2, float scale, float scale2, int32_t scaleMode) const +{ + /// Transformation real y,z -> measured y,z + + GPUCA_RTC_SPECIAL_CODE(ref2 = nullptr; scale2 = 0.f;); + + float dy = 0; + float dz = 0; + + if ((scale >= 0.f) || (scaleMode == 1) || (scaleMode == 2)) { + const auto corrYZ = getCorrectionYZatRealYZ(sector, row, realY, realZ); + dy = corrYZ[0]; + dz = corrYZ[1]; + + if (ref) { // scaling was requested + if (scaleMode == 0 && scale > 0.f) { + const auto val = ref->getCorrectionYZatRealYZ(sector, row, realY, realZ); + dy = (dy - val[0]) * scale + val[0]; + dz = (dz - val[1]) * scale + val[1]; + } else if ((scale != 0) && ((scaleMode == 1) || (scaleMode == 2))) { + const auto val = ref->getCorrectionYZatRealYZ(sector, row, realY, realZ); + dy = val[0] * scale + dy; + dz = val[1] * scale + dz; + } + if (ref2 && (scale2 != 0)) { + const auto val = ref2->getCorrectionYZatRealYZ(sector, row, realY, realZ); + dy = val[0] * scale2 + dy; + dz = val[1] * scale2 + dz; + } + } + } + + measuredY = realY - dy; + measuredZ = realZ - dz; + + GPUCA_DEBUG_STREAMER_CHECK(if (o2::utils::DebugStreamer::checkStream(o2::utils::StreamFlags::streamFastTransform)) { + o2::utils::DebugStreamer::instance()->getStreamer("debug_fasttransform", "UPDATE") << o2::utils::DebugStreamer::instance()->getUniqueTreeName("tree_InverseTransformYZtoNominalYZ").data() + << "sector=" << sector + << "row=" << row + << "scale=" << scale + << "real y=" << realY + << "real z=" << realZ + << "measured y=" << measuredY + << "measured z=" << measuredZ + << "\n"; + }) +} + +GPUdi() void TPCFastTransformPOD::InverseTransformYZtoNominalYZ_new(int32_t sector, int32_t row, float realY, float realZ, float& measuredY, float& measuredZ) const +{ + /// Transformation real y,z -> measured y,z + const auto corrYZ = getCorrectionYZatRealYZ(sector, row, realY, realZ); + measuredY = realY - corrYZ[0]; + measuredZ = realZ - corrYZ[1]; + + GPUCA_DEBUG_STREAMER_CHECK(if (o2::utils::DebugStreamer::checkStream(o2::utils::StreamFlags::streamFastTransform)) { + o2::utils::DebugStreamer::instance()->getStreamer("debug_fasttransform", "UPDATE") << o2::utils::DebugStreamer::instance()->getUniqueTreeName("tree_InverseTransformYZtoNominalYZ").data() + << "sector=" << sector + << "row=" << row + << "real y=" << realY + << "real z=" << realZ + << "measured y=" << measuredY + << "measured z=" << measuredZ + << "\n"; + }) +} + +GPUdi() void TPCFastTransformPOD::InverseTransformXYZtoNominalXYZ(int32_t sector, int32_t row, float x, float y, float z, float& nx, float& ny, float& nz, const TPCFastTransformPOD* ref, const TPCFastTransformPOD* ref2, float scale, float scale2, int32_t scaleMode) const +{ + /// Inverse transformation: Transformed X, Y and Z -> X, Y and Z, transformed w/o space charge correction + int32_t row2 = row + 1; + if (row2 >= getGeometry().getNumberOfRows()) { + row2 = row - 1; + } + float nx1, ny1, nz1; // nominal coordinates for row + float nx2, ny2, nz2; // nominal coordinates for row2 + nx1 = getGeometry().getRowInfo(row).x; + nx2 = getGeometry().getRowInfo(row2).x; + InverseTransformYZtoNominalYZ(sector, row, y, z, ny1, nz1, ref, ref2, scale, scale2, scaleMode); + InverseTransformYZtoNominalYZ(sector, row2, y, z, ny2, nz2, ref, ref2, scale, scale2, scaleMode); + float c1 = (nx2 - nx) / (nx2 - nx1); + float c2 = (nx - nx1) / (nx2 - nx1); + nx = x; + ny = (ny1 * c1 + ny2 * c2); + nz = (nz1 * c1 + nz2 * c2); +} + +GPUdi() void TPCFastTransformPOD::InverseTransformXYZtoNominalXYZ_new(int32_t sector, int32_t row, float x, float y, float z, float& nx, float& ny, float& nz) const +{ + /// Inverse transformation: Transformed X, Y and Z -> X, Y and Z, transformed w/o space charge correction + int32_t row2 = row + 1; + if (row2 >= getGeometry().getNumberOfRows()) { + row2 = row - 1; + } + float nx1, ny1, nz1; // nominal coordinates for row + float nx2, ny2, nz2; // nominal coordinates for row2 + nx1 = getGeometry().getRowInfo(row).x; + nx2 = getGeometry().getRowInfo(row2).x; + InverseTransformYZtoNominalYZ_new(sector, row, y, z, ny1, nz1); + InverseTransformYZtoNominalYZ_new(sector, row2, y, z, ny2, nz2); + float c1 = (nx2 - nx) / (nx2 - nx1); + float c2 = (nx - nx1) / (nx2 - nx1); + nx = x; + ny = (ny1 * c1 + ny2 * c2); + nz = (nz1 * c1 + nz2 * c2); +} + +} // namespace gpu +} // namespace o2 + +#endif diff --git a/GPU/TPCFastTransformation/TPCFastTransformationLinkDef_O2.h b/GPU/TPCFastTransformation/TPCFastTransformationLinkDef_O2.h index f1872549a46aa..0247bbbfbb65b 100644 --- a/GPU/TPCFastTransformation/TPCFastTransformationLinkDef_O2.h +++ b/GPU/TPCFastTransformation/TPCFastTransformationLinkDef_O2.h @@ -93,5 +93,6 @@ #pragma link C++ struct o2::gpu::MultivariatePolynomialContainer + ; #pragma link C++ struct o2::gpu::NDPiecewisePolynomialContainer + ; #pragma link C++ struct o2::gpu::TPCSlowSpaceChargeCorrection + ; +#pragma link C++ class o2::gpu::TPCFastTransformPOD + ; #endif