diff --git a/GPU/TPCFastTransformation/CMakeLists.txt b/GPU/TPCFastTransformation/CMakeLists.txt index b338e1492cc6c..c1c155805117e 100644 --- a/GPU/TPCFastTransformation/CMakeLists.txt +++ b/GPU/TPCFastTransformation/CMakeLists.txt @@ -25,6 +25,7 @@ set(SRCS TPCFastSpaceChargeCorrection.cxx TPCFastSpaceChargeCorrectionMap.cxx TPCFastTransform.cxx + TPCFastTransformPOD.cxx CorrectionMapsHelper.cxx ) diff --git a/GPU/TPCFastTransformation/TPCFastSpaceChargeCorrection.h b/GPU/TPCFastTransformation/TPCFastSpaceChargeCorrection.h index ecbc3c4ad73f5..5560b7d8d8fef 100644 --- a/GPU/TPCFastTransformation/TPCFastSpaceChargeCorrection.h +++ b/GPU/TPCFastTransformation/TPCFastSpaceChargeCorrection.h @@ -37,6 +37,8 @@ namespace gpu /// class TPCFastSpaceChargeCorrection : public FlatObject { + friend class TPCFastTransformPOD; + public: /// /// \brief The struct contains necessary info for TPC padrow diff --git a/GPU/TPCFastTransformation/TPCFastTransform.cxx b/GPU/TPCFastTransformation/TPCFastTransform.cxx index 6424292a6035a..855c63cc2fe2e 100644 --- a/GPU/TPCFastTransformation/TPCFastTransform.cxx +++ b/GPU/TPCFastTransformation/TPCFastTransform.cxx @@ -37,7 +37,7 @@ using namespace GPUCA_NAMESPACE::gpu; TPCFastTransform::TPCFastTransform() - : FlatObject(), mTimeStamp(0), mCorrection(), mApplyCorrection(1), mT0(0.f), mVdrift(0.f), mVdriftCorrY(0.f), mLdriftCorr(0.f), mTOFcorr(0.f), mPrimVtxZ(0.f), mLumi(0.f), mLumiError(0.f), mLumiScaleFactor(1.0f) + : FlatObject(), mTimeStamp(0), mCorrection(), mApplyCorrection(1), mT0(0.f), mVdrift(0.f), mVdriftCorrY(0.f), mLdriftCorr(0.f), mTOFcorr(0.f), mPrimVtxZ(0.f), mLumiIDC(0.f), mLumi(0.f), mLumiError(0.f), mLumiScaleFactor(1.0f) { // Default Constructor: creates an empty uninitialized object } @@ -58,6 +58,7 @@ void TPCFastTransform::cloneFromObject(const TPCFastTransform& obj, char* newFla mLdriftCorr = obj.mLdriftCorr; mTOFcorr = obj.mTOFcorr; mPrimVtxZ = obj.mPrimVtxZ; + mLumiIDC = obj.mLumiIDC; mLumi = obj.mLumi; mLumiError = obj.mLumiError; mLumiScaleFactor = obj.mLumiScaleFactor; @@ -108,6 +109,7 @@ void TPCFastTransform::startConstruction(const TPCFastSpaceChargeCorrection& cor mLdriftCorr = 0.f; mTOFcorr = 0.f; mPrimVtxZ = 0.f; + mLumiIDC = 0.f; mLumi = 0.f; mLumiError = 0.f; mLumiScaleFactor = 1.f; @@ -158,6 +160,7 @@ void TPCFastTransform::print() const LOG(info) << "mLdriftCorr = " << mLdriftCorr; LOG(info) << "mTOFcorr = " << mTOFcorr; LOG(info) << "mPrimVtxZ = " << mPrimVtxZ; + LOG(info) << "mLumiIDC = " << mLumiIDC; LOG(info) << "mLumi = " << mLumi; LOG(info) << "mLumiError = " << mLumiError; LOG(info) << "mLumiScaleFactor = " << mLumiScaleFactor; diff --git a/GPU/TPCFastTransformation/TPCFastTransform.h b/GPU/TPCFastTransformation/TPCFastTransform.h index fe52bffe14acc..c431f8459367c 100644 --- a/GPU/TPCFastTransformation/TPCFastTransform.h +++ b/GPU/TPCFastTransformation/TPCFastTransform.h @@ -162,6 +162,7 @@ class TPCFastTransform : public FlatObject /// Set Lumi info void setLumi(float l) { mLumi = l; } + void setLumiIDC(float l) { mLumiIDC = l; } void setLumiError(float e) { mLumiError = e; } void setLumiScaleFactor(float s) { mLumiScaleFactor = s; } @@ -250,9 +251,15 @@ class TPCFastTransform : public FlatObject /// Return TOF correction (vdrift / C) GPUd() float getTOFCorr() const { return mLdriftCorr; } + /// Return nominal PV Z position correction (vdrift / C) + GPUd() float getPrimVtxZ() const { return mPrimVtxZ; } + /// Return map lumi GPUd() float getLumi() const { return mLumi; } + /// Return map lumiIDC + GPUd() float getLumiIDC() const { return mLumiIDC; } + /// Return map lumi error GPUd() float getLumiError() const { return mLumiError; } @@ -332,6 +339,7 @@ class TPCFastTransform : public FlatObject float mPrimVtxZ; ///< Z of the primary vertex, needed for the Time-Of-Flight correction + float mLumiIDC; ///< luminosity estimator in IDC units float mLumi; ///< luminosity estimator float mLumiError; ///< error on luminosity float mLumiScaleFactor; ///< user correction factor for lumi (e.g. normalization, efficiency correction etc.) @@ -342,7 +350,7 @@ class TPCFastTransform : public FlatObject GPUd() void TransformInternal(int slice, int row, float& u, float& v, float& x, const TPCFastTransform* ref, const TPCFastTransform* ref2, float scale, float scale2, int scaleMode) const; #ifndef GPUCA_ALIROOT_LIB - ClassDefNV(TPCFastTransform, 3); + ClassDefNV(TPCFastTransform, 4); #endif }; diff --git a/GPU/TPCFastTransformation/TPCFastTransformPOD.cxx b/GPU/TPCFastTransformation/TPCFastTransformPOD.cxx new file mode 100644 index 0000000000000..2435bcfb52e0d --- /dev/null +++ b/GPU/TPCFastTransformation/TPCFastTransformPOD.cxx @@ -0,0 +1,243 @@ +// 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" +#include "GPUDebugStreamer.h" +#if !defined(GPUCA_GPUCODE) +#include +#endif + +namespace GPUCA_NAMESPACE +{ +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 slice = 0; slice < origCorr.mGeo.getNumberOfSlices(); slice++) { + for (int row = 0; row < NROWS; row++) { + const auto& spline = origCorr.getSpline(slice, 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.mApplyCorrections = 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 row = 0; row < NROWS; row++) { + podMap.mRowInfo[row] = origCorr.getRowInfo(row); // dataOffsetBytes will be modified later + } + for (int slice = 0; slice < TPCFastTransformGeo::getNumberOfSlices(); slice++) { + podMap.mSliceInfo[slice] = origCorr.getSliceInfo(slice); + for (int row = 0; row < NROWS; row++) { + podMap.mSliceRowInfo[NROWS * slice + row] = origCorr.getSliceRowInfo(slice, row); + } + } + podMap.mInterpolationSafetyMargin = origCorr.fInterpolationSafetyMargin; + podMap.mTimeStamp = origCorr.mTimeStamp; + // + // init data members coming from the TPCFastTrasform + podMap.mVdrift = 0.; + podMap.mT0 = 0.; + podMap.mVdriftCorrY = 0.; + podMap.mLdriftCorr = 0.; + podMap.mTOFcorr = 0.; + podMap.mPrimVtxZ = 0.; + // copy fixed size data --- end + + size_t nextDynOffs = alignOffset(sizeof(TPCFastTransformPOD)); + // copy slice 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 slice = 0; slice < origCorr.mGeo.getNumberOfSlices(); slice++) { + podMap.mSplineDataOffsets[slice][is] = nextDynOffs; + size_t rowDataOffs = 0; + for (int row = 0; row < NROWS; row++) { + const auto& spline = origCorr.getSpline(slice, row); + const float* dataOr = origCorr.getSplineData(slice, 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 slice:{} row:{} to offset {}", nPar, is, slice, 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 slice{}/row{} to {}, overflowing the buffer of size {}", nbcopy, is, slice, row, nextDynOffs, buffSize)); + } + std::memcpy(data, dataOr, nbcopy); + podMap.getRowInfo(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 podMap; +} + +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(); + podMap.mVdriftCorrY = src.getVdriftCorrY(); + podMap.mLdriftCorr = src.getLdriftCorr(); + podMap.mTOFcorr = src.getTOFCorr(); + podMap.mPrimVtxZ = src.getPrimVtxZ(); + // copy fixed size data --- end + return podMap; +} + +bool TPCFastTransformPOD::test(const TPCFastSpaceChargeCorrection& origCorr, int npoints) const +{ + if (npoints < 1) { + return false; + } + std::vector slice, row; + std::vector u, v, dxO, duO, dvO, dxP, duP, dvP, corrXO, corrXP, nomUO, nomVO, nomUP, nomVP; + slice.reserve(npoints); + row.reserve(npoints); + u.reserve(npoints); + v.reserve(npoints); + dxO.resize(npoints); + duO.resize(npoints); + dvO.resize(npoints); + corrXO.resize(npoints); + nomUO.resize(npoints); + nomVO.resize(npoints); + dxP.resize(npoints); + duP.resize(npoints); + dvP.resize(npoints); + corrXP.resize(npoints); + nomUP.resize(npoints); + nomVP.resize(npoints); + + for (int i = 0; i < npoints; i++) { + slice.push_back(gRandom->Integer(NSLICES)); + row.push_back(gRandom->Integer(NROWS)); + u.push_back(gRandom->Rndm() * 15); + v.push_back(gRandom->Rndm() * 200); + } + 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++) { + origCorr.getCorrection(slice[i], row[i], u[i], v[i], dxO[i], duO[i], dvO[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++) { + origCorr.getCorrectionInvCorrectedX(slice[i], row[i], u[i], v[i], corrXO[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++) { + origCorr.getCorrectionInvUV(slice[i], row[i], u[i], v[i], nomUO[i], nomVO[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++) { + this->getCorrection(slice[i], row[i], u[i], v[i], dxP[i], duP[i], dvP[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++) { + this->getCorrectionInvCorrectedX(slice[i], row[i], u[i], v[i], corrXP[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++) { + this->getCorrectionInvUV(slice[i], row[i], u[i], v[i], nomUP[i], nomVP[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 (dxO[i] != dxP[i] || duO[i] != duP[i] || dvO[i] != dvP[i]) { + ndiff[0]++; + } + if (corrXO[i] != corrXP[i]) { + ndiff[1]++; + } + if (nomUO[i] != nomUP[i] || nomVO[i] != nomVP[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 GPUCA_NAMESPACE diff --git a/GPU/TPCFastTransformation/TPCFastTransformPOD.h b/GPU/TPCFastTransformation/TPCFastTransformPOD.h new file mode 100644 index 0000000000000..dd7e42b000c5a --- /dev/null +++ b/GPU/TPCFastTransformation/TPCFastTransformPOD.h @@ -0,0 +1,807 @@ +// 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 the 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 GPUCA_NAMESPACE +{ +namespace gpu +{ + +class TPCFastTransformPOD +{ + public: + using SplineType = TPCFastSpaceChargeCorrection::SplineType; + using RowInfo = TPCFastSpaceChargeCorrection::RowInfo; + using RowActiveArea = TPCFastSpaceChargeCorrection::RowActiveArea; + using SliceRowInfo = TPCFastSpaceChargeCorrection::SliceRowInfo; + using SliceInfo = TPCFastSpaceChargeCorrection::SliceInfo; + + /// convert prefilled buffer to TPCFastTransformPOD + GPUd() static const TPCFastTransformPOD& get(const char* head) { return *reinterpret_cast(head); } + + GPUd() int getApplyCorrections() const { return mApplyCorrections; } + GPUd() void setApplyCorrections(bool v = true) { mApplyCorrections = v; } + + /// _______________ high level methods a la TPCFastTransform _______________________ + /// + GPUd() void Transform(int slice, int row, float pad, float time, float& x, float& y, float& z, float vertexTime = 0) const; + GPUd() void TransformXYZ(int slice, int row, float& x, float& y, float& z) const; + + /// Transformation in the time frame + GPUd() void TransformInTimeFrame(int slice, int row, float pad, float time, float& x, float& y, float& z, float maxTimeBin) const; + GPUd() void TransformInTimeFrame(int slice, float time, float& z, float maxTimeBin) const; + + /// Inverse transformation + GPUd() void InverseTransformInTimeFrame(int slice, int row, float /*x*/, float y, float z, float& pad, float& time, float maxTimeBin) const; + + /// Inverse transformation: Transformed Y and Z -> transformed X + GPUd() void InverseTransformYZtoX(int slice, int 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(int slice, int 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(int slice, int row, float x, float y, float z, float& nx, float& ny, float& nz) const; + + /// Ideal transformation with Vdrift only - without calibration + GPUd() void TransformIdeal(int slice, int row, float pad, float time, float& x, float& y, float& z, float vertexTime) const; + GPUd() void TransformIdealZ(int slice, float time, float& z, float vertexTime) const; + + GPUd() void convPadTimeToUV(int slice, int row, float pad, float time, float& u, float& v, float vertexTime) const; + GPUd() void convPadTimeToUVinTimeFrame(int slice, int row, float pad, float time, float& u, float& v, float maxTimeBin) const; + GPUd() void convTimeToVinTimeFrame(int slice, float time, float& v, float maxTimeBin) const; + + GPUd() void convUVtoPadTime(int slice, int row, float u, float v, float& pad, float& time, float vertexTime) const; + GPUd() void convUVtoPadTimeInTimeFrame(int slice, int row, float u, float v, float& pad, float& time, float maxTimeBin) const; + GPUd() void convVtoTime(float v, float& time, float vertexTime) const; + + GPUd() float convTimeToZinTimeFrame(int slice, float time, float maxTimeBin) const; + GPUd() float convZtoTimeInTimeFrame(int slice, float z, float maxTimeBin) const; + GPUd() float convDeltaTimeToDeltaZinTimeFrame(int slice, float deltaTime) const; + GPUd() float convDeltaZtoDeltaTimeInTimeFrame(int slice, float deltaZ) const; + GPUd() float convDeltaZtoDeltaTimeInTimeFrameAbs(float deltaZ) const; + GPUd() float convZOffsetToVertexTime(int slice, float zOffset, float maxTimeBin) const; + GPUd() float convVertexTimeToZOffset(int slice, float vertexTime, float maxTimeBin) const; + + GPUd() void getTOFcorrection(int slice, int row, float x, float y, float z, float& dz) const; + + /// _______________ methods a la TPCFastSpaceChargeCorrection: cluster correction _______________________ + /// + GPUd() int getCorrection(int slice, int row, float u, float v, float& dx, float& du, float& dv) const; + + /// temporary method with the an way of calculating 2D spline + GPUd() int getCorrectionOld(int slice, int row, float u, float v, float& dx, float& du, float& dv) const; + + /// inverse correction: Corrected U and V -> coorrected X + GPUd() void getCorrectionInvCorrectedX(int slice, int row, float corrU, float corrV, float& corrX) const; + + /// inverse correction: Corrected U and V -> uncorrected U and V + GPUd() void getCorrectionInvUV(int slice, int row, float corrU, float corrV, float& nomU, float& nomV) const; + + /// maximal possible drift length of the active area + GPUd() float getMaxDriftLength(int slice, int row, float pad) const; + + /// maximal possible drift length of the active area + GPUd() float getMaxDriftLength(int slice, int row) const { return getSliceRowInfo(slice, row).activeArea.vMax; } + + /// maximal possible drift length of the active area + GPUd() float getMaxDriftLength(int slice) const { return getSliceInfo(slice).vMax; } + + /// maximal possible drift time of the active area + GPUd() float getMaxDriftTime(int slice, int row, float pad) const; + + /// maximal possible drift time of the active area + GPUd() float getMaxDriftTime(int slice, int row) const; + + /// maximal possible drift time of the active area + GPUd() float getMaxDriftTime(int slice) const; + + /// _______________ Utilities _______________________________________________ + + /// shrink u,v coordinats to the TPC row area +/- fkInterpolationSafetyMargin + GPUd() void schrinkUV(int slice, int row, float& u, float& v) const; + + /// shrink corrected u,v coordinats to the TPC row area +/- fkInterpolationSafetyMargin + GPUd() void schrinkCorrectedUV(int slice, int row, float& corrU, float& corrV) const; + + /// convert u,v to internal grid coordinates + GPUd() void convUVtoGrid(int slice, int row, float u, float v, float& gridU, float& gridV) const; + + /// convert u,v to internal grid coordinates + GPUd() void convGridToUV(int slice, int row, float gridU, float gridV, float& u, float& v) const; + + /// convert corrected u,v to internal grid coordinates + GPUd() void convCorrectedUVtoGrid(int slice, int row, float cu, float cv, float& gridU, float& gridV) const; + + /// TPC geometry information + GPUd() const TPCFastTransformGeo& getGeometry() const { return mGeo; } + + /// 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; } + + /// Gives the interpolation safety marging around the TPC row. + GPUd() float getInterpolationSafetyMargin() const { return mInterpolationSafetyMargin; } + + /// 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 VdriftCorrY in time_bin / cn + GPUd() float getVdriftCorrY() const { return mVdriftCorrY; } + + /// Return LdriftCorr offset in cm + GPUd() float getLdriftCorr() const { return mLdriftCorr; } + + /// Return TOF correction + GPUd() float getTOFCorr() const { return mTOFcorr; } + + /// Return nominal PV Z position + GPUd() float getPrimVtxZ() const { return mPrimVtxZ; } + + /// Sets the time stamp of the current calibaration + GPUd() void setTimeStamp(long int v) { mTimeStamp = v; } + + /// Set safety marging for the interpolation around the TPC row. + /// Outside of this area the interpolation returns the boundary values. + GPUd() void setInterpolationSafetyMargin(float val) { mInterpolationSafetyMargin = val; } + + /// Gives TPC row info + GPUd() const RowInfo& getRowInfo(int row) const { return mRowInfo[row]; } + + /// Gives TPC slice info + GPUd() const SliceInfo& getSliceInfo(int slice) const { return mSliceInfo[slice]; } + + /// Gives a reference to a spline + GPUd() const SplineType& getSpline(int slice, int row) const { return *reinterpret_cast(getThis() + getScenarioOffset(getRowInfo(row).splineScenarioID)); } + + /// Gives pointer to spline data + GPUd() const float* getSplineData(int slice, int row, int iSpline = 0) const { return reinterpret_cast(getThis() + mSplineDataOffsets[slice][iSpline] + getRowInfo(row).dataOffsetBytes[iSpline]); } + + /// Gives TPC slice & row info + GPUd() const SliceRowInfo& getSliceRowInfo(int slice, int row) const { return mSliceRowInfo[NROWS * slice + row]; } + +#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, int npoints = 100000) const { return test(src.getCorrection(), npoints); } + bool test(const TPCFastSpaceChargeCorrection& origCorr, int npoints = 100000) const; +#endif + + static constexpr int NROWS = 152; + static constexpr int NSLICES = TPCFastTransformGeo::getNumberOfSlices(); + static constexpr int NSplineIDs = 3; ///< number of spline data sets for each slice/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 + + /// Gives non-const TPC row info + GPUd() RowInfo& getRowInfo(int row) { return mRowInfo[row]; } + + ///< 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]; } + + GPUd() void TransformInternal(int slice, int row, float& u, float& v, float& x) const; + + bool mApplyCorrections{}; ///< 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::getNumberOfSlices()][NSplineIDs]; ///< start of data for each slice and iSpline data + long int mTimeStamp{}; ///< time stamp of the current calibration + float mInterpolationSafetyMargin{0.1f}; // 10% area around the TPC row. Outside of this area the interpolation returns the boundary values. + + /// _____ Parameters for drift length calculation ____ + /// + /// t = (float) time bin, y = global y + /// + /// L(t,y) = (t-mT0)*(mVdrift + mVdriftCorrY*y ) + mLdriftCorr ____ + /// + float mT0; ///< T0 in [time bin] + float mVdrift; ///< VDrift in [cm/time bin] + float mVdriftCorrY; ///< VDrift correction for global Y[cm] in [1/time bin] + float mLdriftCorr; ///< drift length correction in [cm] + + /// A coefficient for Time-Of-Flight correction: drift length -= EstimatedDistanceToVtx[cm]*mTOFcorr + /// + /// Since this correction requires a knowledge of the spatial position, it is appied after mCorrection, + /// not on the drift length but directly on V coordinate. + /// + /// mTOFcorr == mVdrift/(speed of light) + /// + float mTOFcorr; + + float mPrimVtxZ; ///< Z of the primary vertex, needed for the Time-Of-Flight correction + + TPCFastTransformGeo mGeo; ///< TPC geometry information + TPCFastSpaceChargeCorrection::SliceInfo mSliceInfo[TPCFastTransformGeo::getNumberOfSlices()]; ///< SliceInfo array + RowInfo mRowInfo[NROWS]; + SliceRowInfo mSliceRowInfo[NROWS * TPCFastTransformGeo::getNumberOfSlices()]; + + ClassDefNV(TPCFastTransformPOD, 0); +}; + +GPUdi() int TPCFastTransformPOD::getCorrection(int slice, int row, float u, float v, float& dx, float& du, float& dv) const +{ + const SplineType& spline = getSpline(slice, row); + const float* splineData = getSplineData(slice, row); + float gridU = 0, gridV = 0; + convUVtoGrid(slice, row, u, v, gridU, gridV); + float dxuv[3]; + spline.interpolateU(splineData, gridU, gridV, dxuv); + dx = dxuv[0]; + du = dxuv[1]; + dv = dxuv[2]; + return 0; +} + +GPUdi() int TPCFastTransformPOD::getCorrectionOld(int slice, int row, float u, float v, float& dx, float& du, float& dv) const +{ + const SplineType& spline = getSpline(slice, row); + const float* splineData = getSplineData(slice, row); + float gridU = 0, gridV = 0; + convUVtoGrid(slice, row, u, v, gridU, gridV); + float dxuv[3]; + spline.interpolateUold(splineData, gridU, gridV, dxuv); + dx = dxuv[0]; + du = dxuv[1]; + dv = dxuv[2]; + return 0; +} + +GPUdi() void TPCFastTransformPOD::getCorrectionInvCorrectedX(int slice, int row, float corrU, float corrV, float& x) const +{ + float gridU, gridV; + convCorrectedUVtoGrid(slice, row, corrU, corrV, gridU, gridV); + + const Spline2D& spline = reinterpret_cast&>(getSpline(slice, row)); + const float* splineData = getSplineData(slice, row, 1); + float dx = 0; + spline.interpolateU(splineData, gridU, gridV, &dx); + x = mGeo.getRowInfo(row).x + dx; +} + +GPUdi() void TPCFastTransformPOD::getCorrectionInvUV(int slice, int row, float corrU, float corrV, float& nomU, float& nomV) const +{ + float gridU, gridV; + convCorrectedUVtoGrid(slice, row, corrU, corrV, gridU, gridV); + + const Spline2D& spline = reinterpret_cast&>(getSpline(slice, row)); + const float* splineData = getSplineData(slice, row, 2); + + float duv[2]; + spline.interpolateU(splineData, gridU, gridV, duv); + nomU = corrU - duv[0]; + nomV = corrV - duv[1]; +} + +GPUdi() float TPCFastTransformPOD::getMaxDriftLength(int slice, int row, float pad) const +{ + const RowActiveArea& area = getSliceRowInfo(slice, row).activeArea; + const float* c = area.maxDriftLengthCheb; + float x = -1.f + 2.f * pad / mGeo.getRowInfo(row).maxPad; + float y = c[0] + c[1] * x; + float f0 = 1.f; + float f1 = x; + x *= 2.f; + for (int i = 2; i < 5; i++) { + double f = x * f1 - f0; + y += c[i] * f; + f0 = f1; + f1 = f; + } + return y; +} + +GPUdi() void TPCFastTransformPOD::schrinkUV(int slice, int row, float& u, float& v) const +{ + /// shrink u,v coordinats to the TPC row area +/- mInterpolationSafetyMargin + + float uWidth05 = mGeo.getRowInfo(row).getUwidth() * (0.5f + mInterpolationSafetyMargin); + float vWidth = mGeo.getTPCzLength(slice); + + if (u < -uWidth05) { + u = -uWidth05; + } + if (u > uWidth05) { + u = uWidth05; + } + if (v < -0.1f * vWidth) { + v = -0.1f * vWidth; + } + if (v > 1.1f * vWidth) { + v = 1.1f * vWidth; + } +} + +GPUdi() void TPCFastTransformPOD::schrinkCorrectedUV(int slice, int row, float& corrU, float& corrV) const +{ + /// shrink corrected u,v coordinats to the TPC row area +/- mInterpolationSafetyMargin + + const SliceRowInfo& sliceRowInfo = getSliceRowInfo(slice, row); + + float uMargin = mInterpolationSafetyMargin * mGeo.getRowInfo(row).getUwidth(); + float vMargin = mInterpolationSafetyMargin * mGeo.getTPCzLength(slice); + + if (corrU < sliceRowInfo.activeArea.cuMin - uMargin) { + corrU = sliceRowInfo.activeArea.cuMin - uMargin; + } + + if (corrU > sliceRowInfo.activeArea.cuMax + uMargin) { + corrU = sliceRowInfo.activeArea.cuMax + uMargin; + } + + if (corrV < 0.f - vMargin) { + corrV = 0.f - vMargin; + } + + if (corrV > sliceRowInfo.activeArea.cvMax + vMargin) { + corrV = sliceRowInfo.activeArea.cvMax + vMargin; + } +} + +GPUdi() void TPCFastTransformPOD::convUVtoGrid(int slice, int row, float u, float v, float& gu, float& gv) const +{ + // TODO optimise !!! + gu = 0.f; + gv = 0.f; + + schrinkUV(slice, row, u, v); + + const SliceRowInfo& info = getSliceRowInfo(slice, row); + const SplineType& spline = getSpline(slice, row); + + float su0 = 0.f, sv0 = 0.f; + mGeo.convUVtoScaledUV(slice, row, u, info.gridV0, su0, sv0); + mGeo.convUVtoScaledUV(slice, row, u, v, gu, gv); + + gv = (gv - sv0) / (1.f - sv0); + gu *= spline.getGridX1().getUmax(); + gv *= spline.getGridX2().getUmax(); +} + +GPUdi() void TPCFastTransformPOD::convGridToUV(int slice, int row, float gridU, float gridV, float& u, float& v) const +{ + // TODO optimise + /// convert u,v to internal grid coordinates + float su0 = 0.f, sv0 = 0.f; + const SliceRowInfo& info = getSliceRowInfo(slice, row); + const SplineType& spline = getSpline(slice, row); + mGeo.convUVtoScaledUV(slice, row, 0.f, info.gridV0, su0, sv0); + float su = gridU / spline.getGridX1().getUmax(); + float sv = sv0 + gridV / spline.getGridX2().getUmax() * (1.f - sv0); + mGeo.convScaledUVtoUV(slice, row, su, sv, u, v); +} + +GPUdi() void TPCFastTransformPOD::convCorrectedUVtoGrid(int slice, int row, float corrU, float corrV, float& gridU, float& gridV) const +{ + schrinkCorrectedUV(slice, row, corrU, corrV); + + const SliceRowInfo& sliceRowInfo = getSliceRowInfo(slice, row); + + gridU = (corrU - sliceRowInfo.gridCorrU0) * sliceRowInfo.scaleCorrUtoGrid; + gridV = (corrV - sliceRowInfo.gridCorrV0) * sliceRowInfo.scaleCorrVtoGrid; +} + +#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(); + return create(destVector.data(), destVector.size(), src); +} + +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::Transform(int slice, int 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 slice + /// taking calibration + alignment into account. + const TPCFastTransformGeo::RowInfo& rowInfo = getGeometry().getRowInfo(row); + x = rowInfo.x; + float u = 0, v = 0; + convPadTimeToUV(slice, row, pad, time, u, v, vertexTime); + TransformInternal(slice, row, u, v, x); + getGeometry().convUVtoLocal(slice, u, v, y, z); + float dzTOF = 0; + getTOFcorrection(slice, row, x, y, z, dzTOF); + z += dzTOF; +} + +GPUdi() void TPCFastTransformPOD::TransformXYZ(int slice, int row, float& x, float& y, float& z) const +{ + float u, v; + getGeometry().convLocalToUV(slice, y, z, u, v); + TransformInternal(slice, row, u, v, x); + getGeometry().convUVtoLocal(slice, u, v, y, z); + float dzTOF = 0; + getTOFcorrection(slice, row, x, y, z, dzTOF); + z += dzTOF; +} + +GPUdi() void TPCFastTransformPOD::TransformInTimeFrame(int slice, float time, float& z, float maxTimeBin) const +{ + float v = 0; + convTimeToVinTimeFrame(slice, time, v, maxTimeBin); + getGeometry().convVtoLocal(slice, v, z); +} + +GPUdi() void TPCFastTransformPOD::TransformInTimeFrame(int slice, int 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; + float u = 0, v = 0; + convPadTimeToUVinTimeFrame(slice, row, pad, time, u, v, maxTimeBin); + getGeometry().convUVtoLocal(slice, u, v, y, z); +} + +GPUdi() void TPCFastTransformPOD::InverseTransformInTimeFrame(int slice, int row, float /*x*/, float y, float z, float& pad, float& time, float maxTimeBin) const +{ + /// Inverse transformation to TransformInTimeFrame + float u = 0, v = 0; + getGeometry().convLocalToUV(slice, y, z, u, v); + convUVtoPadTimeInTimeFrame(slice, row, u, v, pad, time, maxTimeBin); +} + +GPUdi() void TPCFastTransformPOD::InverseTransformYZtoX(int slice, int row, float y, float z, float& x) const +{ + /// Transformation y,z -> x + float u = 0, v = 0; + getGeometry().convLocalToUV(slice, y, z, u, v); + if (mApplyCorrections) { + getCorrectionInvCorrectedX(slice, row, u, v, x); + } else { + x = getGeometry().getRowInfo(row).x; // corrections are disabled + } + 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() + << "slice=" << slice + << "row=" << row + << "y=" << y + << "z=" << z + << "x=" << x + << "v=" << v + << "u=" << u + << "\n"; + }) +} + +GPUdi() void TPCFastTransformPOD::InverseTransformYZtoNominalYZ(int slice, int row, float y, float z, float& ny, float& nz) const +{ + /// Transformation y,z -> x + float u = 0, v = 0, un = 0, vn = 0; + getGeometry().convLocalToUV(slice, y, z, u, v); + if (mApplyCorrections) { + getCorrectionInvUV(slice, row, u, v, un, vn); + } else { + un = u; + vn = v; + } + getGeometry().convUVtoLocal(slice, un, vn, ny, nz); + 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() + << "slice=" << slice + << "row=" << row + << "y=" << y + << "z=" << z + << "ny=" << ny + << "nz=" << nz + << "u=" << u + << "v=" << v + << "un=" << un + << "vn=" << vn + << "\n"; + }) +} + +GPUdi() void TPCFastTransformPOD::InverseTransformXYZtoNominalXYZ(int slice, int 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 + int 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(slice, row, y, z, ny1, nz1); + InverseTransformYZtoNominalYZ(slice, 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); +} + +GPUdi() void TPCFastTransformPOD::TransformInternal(int slice, int row, float& u, float& v, float& x) const +{ + if (mApplyCorrections) { + float dx = 0.f, du = 0.f, dv = 0.f; + getCorrection(slice, row, u, v, dx, du, dv); + GPUCA_DEBUG_STREAMER_CHECK(if (o2::utils::DebugStreamer::checkStream(o2::utils::StreamFlags::streamFastTransform)) { + float ly, lz; + getGeometry().convUVtoLocal(slice, u, v, ly, lz); + float gx, gy, gz; + getGeometry().convLocalToGlobal(slice, x, ly, lz, gx, gy, gz); + float lyT, lzT; + float uCorr = u + du; + float vCorr = v + dv; + float lxT = x + dx; + getGeometry().convUVtoLocal(slice, uCorr, vCorr, lyT, lzT); + float invYZtoX; + InverseTransformYZtoX(slice, row, lyT, lzT, invYZtoX); + float YZtoNominalY; + float YZtoNominalZ; + InverseTransformYZtoNominalYZ(slice, 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 + << "du=" << du + << "dv=" << dv + << "v=" << v + << "u=" << u + << "row=" << row + << "slice=" << slice + // original local coordinates + << "ly=" << ly + << "lz=" << lz + << "lx=" << x + // 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 += dx; + u += du; + v += dv; + } +} + +GPUdi() void TPCFastTransformPOD::TransformIdealZ(int slice, float time, float& z, float vertexTime) const +{ + /// _______________ The main method: cluster transformation _______________________ + /// + /// Transforms time TPC coordinates to local Z withing a slice + /// Ideal transformation: only Vdrift from DCS. + /// No space charge corrections, no time of flight correction + /// + + float v = (time - mT0 - vertexTime) * mVdrift; // drift length cm + getGeometry().convVtoLocal(slice, v, z); +} + +GPUdi() void TPCFastTransformPOD::TransformIdeal(int slice, int 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 slice + /// Ideal transformation: only Vdrift from DCS. + /// No space charge corrections, no time of flight correction + /// + const TPCFastTransformGeo::RowInfo& rowInfo = getGeometry().getRowInfo(row); + x = rowInfo.x; + float u = (pad - 0.5f * rowInfo.maxPad) * rowInfo.padWidth; + float v = (time - mT0 - vertexTime) * mVdrift; // drift length cm + getGeometry().convUVtoLocal(slice, u, v, y, z); +} + +GPUdi() void TPCFastTransformPOD::convPadTimeToUV(int slice, int row, float pad, float time, float& u, float& v, float vertexTime) const +{ + bool sideC = (slice >= getGeometry().getNumberOfSlicesA()); + const TPCFastTransformGeo::RowInfo& rowInfo = getGeometry().getRowInfo(row); + const TPCFastTransformGeo::SliceInfo& sliceInfo = getGeometry().getSliceInfo(slice); + float x = rowInfo.x; + u = (pad - 0.5f * rowInfo.maxPad) * rowInfo.padWidth; + float y = sideC ? -u : u; // pads are mirrorred on C-side + float yLab = y * sliceInfo.cosAlpha + x * sliceInfo.sinAlpha; + v = (time - mT0 - vertexTime) * (mVdrift + mVdriftCorrY * yLab) + mLdriftCorr; // drift length cm +} + +GPUdi() void TPCFastTransformPOD::convTimeToVinTimeFrame(int slice, float time, float& v, float maxTimeBin) const +{ + v = (time - mT0 - maxTimeBin) * mVdrift + mLdriftCorr; // drift length cm + if (slice < getGeometry().getNumberOfSlicesA()) { + v += getGeometry().getTPCzLengthA(); + } else { + v += getGeometry().getTPCzLengthC(); + } +} + +GPUdi() void TPCFastTransformPOD::convPadTimeToUVinTimeFrame(int slice, int row, float pad, float time, float& u, float& v, float maxTimeBin) const +{ + const TPCFastTransformGeo::RowInfo& rowInfo = getGeometry().getRowInfo(row); + u = (pad - 0.5f * rowInfo.maxPad) * rowInfo.padWidth; + convTimeToVinTimeFrame(slice, time, v, maxTimeBin); +} + +GPUdi() float TPCFastTransformPOD::convZOffsetToVertexTime(int slice, float zOffset, float maxTimeBin) const +{ + if (slice < getGeometry().getNumberOfSlicesA()) { + return maxTimeBin - (getGeometry().getTPCzLengthA() + zOffset) / mVdrift; + } else { + return maxTimeBin - (getGeometry().getTPCzLengthC() - zOffset) / mVdrift; + } +} + +GPUdi() float TPCFastTransformPOD::convVertexTimeToZOffset(int slice, float vertexTime, float maxTimeBin) const +{ + if (slice < getGeometry().getNumberOfSlicesA()) { + return (maxTimeBin - vertexTime) * mVdrift - getGeometry().getTPCzLengthA(); + } else { + return -((maxTimeBin - vertexTime) * mVdrift - getGeometry().getTPCzLengthC()); + } +} + +GPUdi() void TPCFastTransformPOD::convUVtoPadTime(int slice, int row, float u, float v, float& pad, float& time, float vertexTime) const +{ + bool sideC = (slice >= getGeometry().getNumberOfSlicesA()); + + const TPCFastTransformGeo::RowInfo& rowInfo = getGeometry().getRowInfo(row); + const TPCFastTransformGeo::SliceInfo& sliceInfo = getGeometry().getSliceInfo(slice); + + pad = u / rowInfo.padWidth + 0.5f * rowInfo.maxPad; + + float x = rowInfo.x; + float y = sideC ? -u : u; // pads are mirrorred on C-side + float yLab = y * sliceInfo.cosAlpha + x * sliceInfo.sinAlpha; + time = mT0 + vertexTime + (v - mLdriftCorr) / (mVdrift + mVdriftCorrY * yLab); +} + +GPUdi() void TPCFastTransformPOD::convVtoTime(float v, float& time, float vertexTime) const +{ + float yLab = 0.f; + time = mT0 + vertexTime + (v - mLdriftCorr) / (mVdrift + mVdriftCorrY * yLab); +} + +GPUdi() void TPCFastTransformPOD::convUVtoPadTimeInTimeFrame(int slice, int row, float u, float v, float& pad, float& time, float maxTimeBin) const +{ + if (slice < getGeometry().getNumberOfSlicesA()) { + v -= getGeometry().getTPCzLengthA(); + } else { + v -= getGeometry().getTPCzLengthC(); + } + const TPCFastTransformGeo::RowInfo& rowInfo = getGeometry().getRowInfo(row); + pad = u / rowInfo.padWidth + 0.5f * rowInfo.maxPad; + time = mT0 + maxTimeBin + (v - mLdriftCorr) / mVdrift; +} + +GPUdi() void TPCFastTransformPOD::getTOFcorrection(int slice, int /*row*/, float x, float y, float z, float& dz) const +{ + // calculate time of flight correction for z coordinate + + bool sideC = (slice >= getGeometry().getNumberOfSlicesA()); + float distZ = z - mPrimVtxZ; + float dv = -GPUCommonMath::Sqrt(x * x + y * y + distZ * distZ) * mTOFcorr; + dz = sideC ? dv : -dv; +} + +GPUdi() float TPCFastTransformPOD::getMaxDriftTime(int slice, int row, float pad) const +{ + /// maximal possible drift time of the active area + float maxL = getMaxDriftLength(slice, row, pad); + bool sideC = (slice >= getGeometry().getNumberOfSlicesA()); + const TPCFastTransformGeo::RowInfo& rowInfo = getGeometry().getRowInfo(row); + const TPCFastTransformGeo::SliceInfo& sliceInfo = getGeometry().getSliceInfo(slice); + float x = rowInfo.x; + float u = (pad - 0.5f * rowInfo.maxPad) * rowInfo.padWidth; + + float y = sideC ? -u : u; // pads are mirrorred on C-side + float yLab = y * sliceInfo.cosAlpha + x * sliceInfo.sinAlpha; + return mT0 + (maxL - mLdriftCorr) / (mVdrift + mVdriftCorrY * yLab); +} + +GPUdi() float TPCFastTransformPOD::getMaxDriftTime(int slice, int row) const +{ + /// maximal possible drift time of the active area + float maxL = getMaxDriftLength(slice, row); + float maxTime = 0.f; + convVtoTime(maxL, maxTime, 0.f); + return maxTime; +} + +GPUdi() float TPCFastTransformPOD::getMaxDriftTime(int slice) const +{ + /// maximal possible drift time of the active area + float maxL = getMaxDriftLength(slice); + float maxTime = 0.f; + convVtoTime(maxL, maxTime, 0.f); + return maxTime; +} + +} // namespace gpu +} // namespace GPUCA_NAMESPACE + +#endif // ALICEO2_GPU_TPCFastTransformPOD_H diff --git a/GPU/TPCFastTransformation/TPCFastTransformationLinkDef_O2.h b/GPU/TPCFastTransformation/TPCFastTransformationLinkDef_O2.h index 4421d44aab0c8..7482cbf52768a 100644 --- a/GPU/TPCFastTransformation/TPCFastTransformationLinkDef_O2.h +++ b/GPU/TPCFastTransformation/TPCFastTransformationLinkDef_O2.h @@ -63,6 +63,7 @@ #pragma link C++ class o2::gpu::TPCFastTransformGeo::RowInfo + ; #pragma link C++ class o2::gpu::TPCFastTransform + ; +#pragma link C++ class o2::gpu::TPCFastTransformPOD; #pragma link C++ class o2::gpu::TPCFastSpaceChargeCorrectionMap + ; #pragma link C++ class o2::gpu::TPCFastSpaceChargeCorrection::RowInfo + ;