diff --git a/Detectors/TPC/calibration/include/TPCCalibration/TPCScaler.h b/Detectors/TPC/calibration/include/TPCCalibration/TPCScaler.h index 03a55a4b732b4..55809b55c514e 100644 --- a/Detectors/TPC/calibration/include/TPCCalibration/TPCScaler.h +++ b/Detectors/TPC/calibration/include/TPCCalibration/TPCScaler.h @@ -28,11 +28,18 @@ Class for storing the scalers which are used to calculate an estimate for the me */ struct TPCScalerWeights { - float getWeight(float deltaTime) const; - bool isValid() const { return (mSamplingTimeMS > 0 && !mWeights.empty()); } - float mSamplingTimeMS = -1; ///< sampling of the stored weights - float mFirstTimeStampMS = 0; ///< first timestamp - std::vector mWeights; ///< stored weights + float getWeight(float deltaTime, int idxZ, Side side) const; + const auto& getWeights(Side side) const { return (side == Side::A) ? mWeightsA : mWeightsC; } + bool isValid(Side side) const { return (mSamplingTimeMS > 0 && !getWeights(side).empty()); } + int getNZPoints() const { return mZPos.size(); } + int getNTimePoints() const { return mWeightsA.size() / mZPos.size(); } + int getIdx(int idxTime, int idxZ) const { return idxTime + idxZ * getNTimePoints(); } + float getZPos(int idxZ) const { return mZPos[idxZ]; } + float mSamplingTimeMS = -1; ///< sampling of the timestamps + float mFirstTimeStampMS = 0; ///< first timestamp + std::vector mWeightsA; ///< stored weights for A-side + std::vector mWeightsC; ///< stored weights for C-side + std::vector mZPos; ///< z positions ClassDefNV(TPCScalerWeights, 1); }; @@ -119,9 +126,12 @@ class TPCScaler /// \return returns numbers of scalers for one ion drift time int getNValuesIonDriftTime() const { return mIonDriftTimeMS / mIntegrationTimeMS + 0.5; } - /// \return returns mean scaler value for last ion drift time + /// \return returns mean scaler value for last ion drift time (in case z dependent weights are used return vector of weights) /// \param timestamp timestamp for which the last values are used to calculate the mean - float getMeanScaler(double timestamp, o2::tpc::Side side) const; + std::vector getMeanScaler(double timestamp, o2::tpc::Side side) const; + + /// \return returns mean over all stored values + float getMeanScaler(o2::tpc::Side side) const; /// \return returns duration in ms for which the scalers are defined float getDurationMS(o2::tpc::Side side) const { return mIntegrationTimeMS * getNValues(side); } @@ -132,24 +142,25 @@ class TPCScaler void clampScalers(float minThreshold, float maxThreshold, Side side); /// setting the weights for the scalers - void setScalerWeights(const TPCScalerWeights& weights, Side side) { (side == o2::tpc::Side::A) ? mScalerWeightsA = weights : mScalerWeightsC = weights; } + void setScalerWeights(const TPCScalerWeights& weights) { mScalerWeights = weights; } /// \return returns stored weights for TPC scalers - const auto& getScalerWeights(o2::tpc::Side side) const { return (side == o2::tpc::Side::A) ? mScalerWeightsA : mScalerWeightsC; } + const auto& getScalerWeights() const { return mScalerWeights; } void useWeights(bool useWeights) { mUseWeights = useWeights; } + std::vector getZPosWeights() const { return mScalerWeights.mZPos; } + private: - float mIonDriftTimeMS{200}; ///< ion drift time in ms - int mRun{}; ///< run for which this object is valid - unsigned int mFirstTFOrbit{}; ///< first TF orbit of the stored scalers - double mTimeStampMS{}; ///< time stamp of the first stored values - float mIntegrationTimeMS{1.}; ///< integration time for each stored value in ms - std::vector mScalerA{}; ///< TPC scaler for A-side - std::vector mScalerC{}; ///< TPC scaler for C-side - TPCScalerWeights mScalerWeightsA{}; ///< weights for the scalers for A-side - TPCScalerWeights mScalerWeightsC{}; ///< weights for the scalers for C-side - bool mUseWeights{false}; ///< use weights when calculating the mean scaler + float mIonDriftTimeMS{200}; ///< ion drift time in ms + int mRun{}; ///< run for which this object is valid + unsigned int mFirstTFOrbit{}; ///< first TF orbit of the stored scalers + double mTimeStampMS{}; ///< time stamp of the first stored values + float mIntegrationTimeMS{1.}; ///< integration time for each stored value in ms + std::vector mScalerA{}; ///< TPC scaler for A-side + std::vector mScalerC{}; ///< TPC scaler for C-side + TPCScalerWeights mScalerWeights{}; ///< weights for the scalers for A-side + bool mUseWeights{false}; ///< use weights when calculating the mean scaler // get index to data for given timestamp int getDataIdx(double timestamp) const { return (timestamp - mTimeStampMS) / mIntegrationTimeMS + 0.5; } diff --git a/Detectors/TPC/calibration/src/CorrectionMapsLoader.cxx b/Detectors/TPC/calibration/src/CorrectionMapsLoader.cxx index d5fa5ad8fd68c..2375fcaf3a0f5 100644 --- a/Detectors/TPC/calibration/src/CorrectionMapsLoader.cxx +++ b/Detectors/TPC/calibration/src/CorrectionMapsLoader.cxx @@ -57,10 +57,19 @@ void CorrectionMapsLoader::extractCCDBInputs(ProcessingContext& pc) } lumiObj = lumiPrev; } - setInstLumi(mInstLumiFactor * (mCTPLumiSource == 0 ? lumiObj.getLumi() : lumiObj.getLumiAlt())); + + const float lumi = mInstLumiFactor * (mCTPLumiSource == 0 ? lumiObj.getLumi() : lumiObj.getLumiAlt()); + const float v[getNBinsLumi()]{lumi, lumi}; + setInstLumi(v); } else if (getLumiScaleType() == 2 && mInstLumiOverride <= 0.) { - float tpcScaler = pc.inputs().get("tpcscaler"); - setInstLumi(mInstLumiFactor * tpcScaler); + std::vector tpcScaler = pc.inputs().get>("tpcscaler"); + + if (tpcScaler.size() != getNBinsLumi()) { + const float v[getNBinsLumi()]{mInstLumiFactor * tpcScaler[0], mInstLumiFactor * tpcScaler[1]}; + setInstLumi(v); + } else { + setInstLumi(tpcScaler.data()); + } } } @@ -170,7 +179,7 @@ bool CorrectionMapsLoader::accountCCDBInputs(const ConcreteDataMatcher& matcher, setMeanLumiRef(mMeanLumiRefOverride); } if (mInstLumiOverride != 0.) { - setInstLumi(mInstLumiOverride, false); + // setInstLumi(mInstLumiOverride, false); } setUpdatedLumi(); int scaleType = getLumiScaleType(); diff --git a/Detectors/TPC/calibration/src/TPCScaler.cxx b/Detectors/TPC/calibration/src/TPCScaler.cxx index d705232bfdd91..3c4240d5cd037 100644 --- a/Detectors/TPC/calibration/src/TPCScaler.cxx +++ b/Detectors/TPC/calibration/src/TPCScaler.cxx @@ -19,6 +19,7 @@ #include #include "Framework/Logger.h" #include "CommonConstants/LHCConstants.h" +#include using namespace o2::tpc; @@ -113,14 +114,14 @@ void TPCScaler::setFromTree(TTree& tpcScalerTree) tpcScalerTree.SetBranchAddress("TPCScaler", nullptr); } -float TPCScaler::getMeanScaler(double timestamp, o2::tpc::Side side) const +std::vector TPCScaler::getMeanScaler(double timestamp, o2::tpc::Side side) const { // index to data buffer const int idxData = getDataIdx(timestamp); const int nVals = getNValuesIonDriftTime(); const int nValues = getNValues(side); if ((nVals == 0) || (nVals > nValues)) { - return -1; + std::vector(-1); LOGP(error, "Empty data provided {}", nValues); } @@ -128,23 +129,38 @@ float TPCScaler::getMeanScaler(double timestamp, o2::tpc::Side side) const const int lastIdx = std::clamp(idxData, nVals, nValues); const int firstIdx = (lastIdx == nValues) ? (nValues - nVals) : std::clamp(idxData - nVals, 0, nValues); - // sump up values from last ion drift time - float sum = 0; - float sumW = 0; - const bool useWeights = mUseWeights && getScalerWeights(side).isValid(); - for (int i = firstIdx; i < lastIdx; ++i) { - float weight = 1; - if (useWeights) { - const double relTSMS = mTimeStampMS + i * mIntegrationTimeMS - timestamp; - weight = getScalerWeights(side).getWeight(relTSMS); + std::vector vWeights; + const bool useWeights = mUseWeights && getScalerWeights().isValid(side); + int izEnd = (useWeights) ? mScalerWeights.getNZPoints() : 1; + vWeights.reserve(izEnd); + for (int iz = 0; iz < izEnd; ++iz) { + // sump up values from last ion drift time + float sum = 0; + float sumW = 0; + for (int i = firstIdx; i < lastIdx; ++i) { + float weight = 1; + // overwrite weight, in case requested + if (useWeights) { + const double relTSMS = mTimeStampMS + i * mIntegrationTimeMS - timestamp; + weight = getScalerWeights().getWeight(relTSMS, iz, side); + } + sum += getScalers(i, side) * weight; + sumW += weight; + } + + if (sumW != 0) { + vWeights.emplace_back(sum / sumW); + } else { + vWeights.emplace_back(0); } - sum += getScalers(i, side) * weight; - sumW += weight; - } - if (sumW != 0) { - return (sum / sumW); } - return 0; + // LOGP(info, "Return weights of size: {} for {} z bins", vWeights, izEnd); + return vWeights; +} + +float TPCScaler::getMeanScaler(o2::tpc::Side side) const +{ + return getScalers(side).empty() ? 0. : (std::accumulate(getScalers(side).begin(), getScalers(side).end(), 0.) / getScalers(side).size()); } void TPCScaler::clampScalers(float minThreshold, float maxThreshold, Side side) @@ -153,12 +169,12 @@ void TPCScaler::clampScalers(float minThreshold, float maxThreshold, Side side) std::transform(std::begin(scaler), std::end(scaler), std::begin(scaler), [minThreshold, maxThreshold](auto val) { return std::clamp(val, minThreshold, maxThreshold); }); } -float TPCScalerWeights::getWeight(float deltaTime) const +float TPCScalerWeights::getWeight(float deltaTime, int idxZ, Side side) const { - const int idx = (deltaTime - mFirstTimeStampMS) / mSamplingTimeMS + 0.5; - if ((idx < 0) || (idx > mWeights.size() - 1)) { + const int idxTime = (deltaTime - mFirstTimeStampMS) / mSamplingTimeMS + 0.5; + if ((idxTime < 0) || (idxTime > getNTimePoints() - 1) || (idxZ < 0) || (idxZ > getNZPoints() - 1)) { // LOGP(info, "Index out of range for deltaTime: {} mFirstTimeStampMS: {} mSamplingTimeMS: {}", deltaTime, mFirstTimeStampMS, mSamplingTimeMS); return 1; } - return mWeights[idx]; + return (side == Side::A) ? mWeightsA[getIdx(idxTime, idxZ)] : mWeightsC[getIdx(idxTime, idxZ)]; } \ No newline at end of file diff --git a/Detectors/TPC/workflow/src/TPCScalerSpec.cxx b/Detectors/TPC/workflow/src/TPCScalerSpec.cxx index fcb16c0fbd119..2ba49ea36156f 100644 --- a/Detectors/TPC/workflow/src/TPCScalerSpec.cxx +++ b/Detectors/TPC/workflow/src/TPCScalerSpec.cxx @@ -62,11 +62,29 @@ class TPCScalerSpec : public Task const auto orbitResetTimeMS = o2::base::GRPGeomHelper::instance().getOrbitResetTimeMS(); const auto firstTFOrbit = pc.services().get().firstTForbit; const double timestamp = orbitResetTimeMS + firstTFOrbit * o2::constants::lhc::LHCOrbitMUS * 0.001; - float scalerA = mTPCScaler.getMeanScaler(timestamp, o2::tpc::Side::A); - float scalerC = mTPCScaler.getMeanScaler(timestamp, o2::tpc::Side::C); - float meanScaler = (scalerA + scalerC) / 2; - LOGP(info, "Publishing TPC scaler: {}", meanScaler); - pc.outputs().snapshot(Output{header::gDataOriginTPC, "TPCSCALER"}, meanScaler); + std::vector scalerA{mTPCScaler.getMeanScaler(timestamp, o2::tpc::Side::A)}; + std::vector scalerC{mTPCScaler.getMeanScaler(timestamp, o2::tpc::Side::C)}; + const std::vector zPosWeights = mTPCScaler.getZPosWeights(); + + if ((scalerA.size() != scalerC.size()) || (scalerA.size() != zPosWeights.size())) { + LOGP(error, "Extracted TPC scalers of different size! Size A-side: {} vs Size C-side: {} vs zPosWeights: {}", scalerA.size(), scalerC.size(), zPosWeights.size()); + } + + const float meanIDC = (mMeanIDCA + mMeanIDCC) / 2; + const float scaleA = mMeanIDCA / meanIDC; + const float scaleC = mMeanIDCC / meanIDC; + LOGP(info, "meanIDC: {} meanIDCA: {} meanIDCC: {} scaleA: {} scaleC: {}", meanIDC, mMeanIDCA, mMeanIDCC, scaleA, scaleC); + std::transform(scalerA.begin(), scalerA.end(), scalerA.begin(), [scaleA](float val) { return val / scaleA; }); + std::transform(scalerC.begin(), scalerC.end(), scalerC.begin(), [scaleC](float val) { return val / scaleC; }); + + // append vector C to A and weights + scalerA.insert(scalerA.end(), scalerC.begin(), scalerC.end()); + scalerA.insert(scalerA.end(), zPosWeights.begin(), zPosWeights.end()); + + for (int i = 0; i < scalerA.size(); ++i) { + LOGP(info, "Publishing TPC scaler: {}: {}", i, scalerA[i]); + } + pc.outputs().snapshot(Output{header::gDataOriginTPC, "TPCSCALER"}, scalerA); } void finaliseCCDB(o2::framework::ConcreteDataMatcher& matcher, void* obj) final @@ -75,22 +93,23 @@ class TPCScalerSpec : public Task if (matcher == ConcreteDataMatcher(o2::header::gDataOriginTPC, "TPCSCALERCCDB", 0)) { LOGP(info, "Updating TPC scaler"); mTPCScaler.setFromTree(*((TTree*)obj)); + mMeanIDCA = mTPCScaler.getMeanScaler(o2::tpc::Side::A); + mMeanIDCC = mTPCScaler.getMeanScaler(o2::tpc::Side::C); + LOGP(info, "Updated mean scaler A {} mean scaler C", mMeanIDCA, mMeanIDCC); if (mIonDriftTimeMS > 0) { LOGP(info, "Setting ion drift time to: {}", mIonDriftTimeMS); mTPCScaler.setIonDriftTimeMS(mIonDriftTimeMS); } - if (mScalerWeights.isValid()) { + if (mScalerWeights.isValid(Side::A)) { LOGP(info, "Setting TPC scaler weights"); - mTPCScaler.setScalerWeights(mScalerWeights, Side::A); - mTPCScaler.setScalerWeights(mScalerWeights, Side::C); + mTPCScaler.setScalerWeights(mScalerWeights); mTPCScaler.useWeights(true); } } if (matcher == ConcreteDataMatcher(o2::header::gDataOriginTPC, "TPCSCALERWCCDB", 0)) { LOGP(info, "Updating TPC scaler weights"); mScalerWeights = *(TPCScalerWeights*)obj; - mTPCScaler.setScalerWeights(mScalerWeights, Side::A); - mTPCScaler.setScalerWeights(mScalerWeights, Side::C); + mTPCScaler.setScalerWeights(mScalerWeights); mTPCScaler.useWeights(true); } } @@ -101,6 +120,8 @@ class TPCScalerSpec : public Task TPCScalerWeights mScalerWeights{}; ///< scaler weights float mIonDriftTimeMS{-1}; ///< ion drift time TPCScaler mTPCScaler; ///< tpc scaler + float mMeanIDCA{0}; + float mMeanIDCC{0}; }; o2::framework::DataProcessorSpec getTPCScalerSpec(bool enableWeights) diff --git a/GPU/TPCFastTransformation/CorrectionMapsHelper.cxx b/GPU/TPCFastTransformation/CorrectionMapsHelper.cxx index cf7fdd8370f78..182285021db1f 100644 --- a/GPU/TPCFastTransformation/CorrectionMapsHelper.cxx +++ b/GPU/TPCFastTransformation/CorrectionMapsHelper.cxx @@ -24,7 +24,9 @@ void CorrectionMapsHelper::clear() mCorrMap = nullptr; mCorrMapRef = nullptr; mUpdatedFlags = 0; - mInstLumi = 0.f; + for (int i = 0; i < getNBinsLumi(); ++i) { + mInstLumi[i] = 0.f; + } mMeanLumi = 0.f; mMeanLumiRef = 0.f; } @@ -78,5 +80,5 @@ void CorrectionMapsHelper::setCorrMapRef(std::unique_ptr&& m) //________________________________________________________ void CorrectionMapsHelper::reportScaling() { - LOGP(info, "Map scaling update: InstLumiOverride={}, LumiScaleType={} -> instLumi={}, meanLumi={} -> LumiScale={}, lumiScaleMode={}", getInstLumiOverride(), getLumiScaleType(), getInstLumi(), getMeanLumi(), getLumiScale(), getLumiScaleMode()); + // LOGP(info, "Map scaling update: InstLumiOverride={}, LumiScaleType={} -> instLumi={}, meanLumi={} -> LumiScale={}, lumiScaleMode={}", getInstLumiOverride(), getLumiScaleType(), getInstLumi(), getMeanLumi(), getLumiScale(), getLumiScaleMode()); } diff --git a/GPU/TPCFastTransformation/CorrectionMapsHelper.h b/GPU/TPCFastTransformation/CorrectionMapsHelper.h index 970f9b4c88a26..6847391ae7a88 100644 --- a/GPU/TPCFastTransformation/CorrectionMapsHelper.h +++ b/GPU/TPCFastTransformation/CorrectionMapsHelper.h @@ -64,10 +64,16 @@ class CorrectionMapsHelper void setCorrMap(GPUCA_NAMESPACE::gpu::TPCFastTransform* m); void setCorrMapRef(GPUCA_NAMESPACE::gpu::TPCFastTransform* m); void reportScaling(); - void setInstLumi(float v, bool report = true) + void setInstLumi(const float v[/* 3 * NLUMIZBINS */], bool report = true) { - if (v != mInstLumi) { - mInstLumi = v; + bool updatelumi = false; + for (int iz = 0; iz < getNBinsLumi(); ++iz) { + if (v[iz] != mInstLumi[iz]) { + mInstLumi[iz] = v[iz]; + updatelumi = true; + } + } + if (updatelumi) { updateLumiScale(report); } } @@ -97,13 +103,27 @@ class CorrectionMapsHelper void updateLumiScale(bool report = true) { - if (mMeanLumi < 0.f || mInstLumi < 0.f) { - mLumiScale = -1.f; + if (mMeanLumi < 0.f || mInstLumi[0] < 0.f) { + for (int i = 0; i < getNBinsLumi(); ++i) { + mLumiScale[i] = -1.f; + } } else if ((mLumiScaleMode == 1) || (mLumiScaleMode == 2)) { - mLumiScale = mMeanLumiRef ? (mInstLumi - mMeanLumi) / mMeanLumiRef : 0.f; - LOGP(debug, "mInstLumi: {} mMeanLumi: {} mMeanLumiRef: {}", mInstLumi, mMeanLumi, mMeanLumiRef); + for (int i = 0; i < NLUMIZBINS; ++i) { + mLumiScale[i] = mMeanLumiRef ? (mInstLumi[i] - mMeanLumi) / mMeanLumiRef : 0.f; + const int idxC = i + NLUMIZBINS; + mLumiScale[idxC] = mMeanLumiRef ? (mInstLumi[idxC] - mMeanLumi) / mMeanLumiRef : 0.f; + const int idxZ = i + 2 * NLUMIZBINS; + mLumiScale[idxZ] = mInstLumi[idxZ]; + LOGP(info, "mInstLumi A: {}, mInstLumi C: {} mMeanLumi: {} mMeanLumiRef: {} for iz: {}", mInstLumi[i], mInstLumi[idxC], mMeanLumi, mMeanLumiRef, mLumiScale[idxZ]); + } } else { - mLumiScale = mMeanLumi ? mInstLumi / mMeanLumi : 0.f; + for (int i = 0; i < NLUMIZBINS; ++i) { + mLumiScale[i] = mMeanLumi ? mInstLumi[i] / mMeanLumi : 0.f; + const int idxC = i + NLUMIZBINS; + mLumiScale[idxC] = mMeanLumi ? mInstLumi[idxC] / mMeanLumi : 0.f; + const int idxZ = i + 2 * NLUMIZBINS; + mLumiScale[idxZ] = mInstLumi[idxZ]; + } } setUpdatedLumi(); if (report) { @@ -111,11 +131,11 @@ class CorrectionMapsHelper } } - GPUd() float getInstLumi() const { return mInstLumi; } + GPUd() const float* getInstLumi() const { return mInstLumi; } GPUd() float getMeanLumi() const { return mMeanLumi; } GPUd() float getMeanLumiRef() const { return mMeanLumiRef; } - GPUd() float getLumiScale() const { return mLumiScale; } + GPUd() const float* getLumiScale() const { return mLumiScale; } GPUd() int getLumiScaleMode() const { return mLumiScaleMode; } bool isUpdated() const { return mUpdatedFlags != 0; } @@ -146,6 +166,10 @@ class CorrectionMapsHelper int getUpdateFlags() const { return mUpdatedFlags; } + static constexpr int getNBinsLumi() { return 3 * NLUMIZBINS; } + + GPUd() static constexpr int getNZBins() { return NLUMIZBINS; } + protected: enum UpdateFlags { MapBit = 0x1, MapRefBit = 0x2, @@ -153,10 +177,11 @@ class CorrectionMapsHelper bool mOwner = false; // is content of pointers owned by the helper int mLumiScaleType = 0; // require CTP Lumi for mInstLumi int mUpdatedFlags = 0; - float mInstLumi = 0.; // instanteneous luminosity (a.u) + constexpr static int NLUMIZBINS = 10; + float mInstLumi[3 * NLUMIZBINS]{0.}; // instanteneous luminosity (a.u) float mMeanLumi = 0.; // mean luminosity of the map (a.u) float mMeanLumiRef = 0.; // mean luminosity of the ref map (a.u) - float mLumiScale = 0.; // precalculated mInstLumi/mMeanLumi + float mLumiScale[3 * NLUMIZBINS]{0.}; // precalculated mInstLumi/mMeanLumi: values, Z-Bins int mLumiScaleMode = 0; // scaling-mode of the correciton maps float mMeanLumiOverride = -1.f; // optional value to override mean lumi float mMeanLumiRefOverride = -1.f; // optional value to override ref mean lumi diff --git a/GPU/TPCFastTransformation/TPCFastTransform.h b/GPU/TPCFastTransformation/TPCFastTransform.h index 1936ae4f19a99..9a2da2659b8ea 100644 --- a/GPU/TPCFastTransformation/TPCFastTransform.h +++ b/GPU/TPCFastTransformation/TPCFastTransform.h @@ -176,8 +176,8 @@ class TPCFastTransform : public FlatObject /// Transforms raw TPC coordinates to local XYZ withing a slice /// taking calibration + alignment into account. /// - GPUd() void Transform(int slice, int row, float pad, float time, float& x, float& y, float& z, float vertexTime = 0, const TPCFastTransform* ref = nullptr, float scale = 0.f, int scaleMode = 0) const; - GPUd() void TransformXYZ(int slice, int row, float& x, float& y, float& z, const TPCFastTransform* ref = nullptr, float scale = 0.f, int scaleMode = 0) const; + GPUd() void Transform(int slice, int row, float pad, float time, float& x, float& y, float& z, float vertexTime = 0, const TPCFastTransform* ref = nullptr, const float* scale = nullptr, int scaleMode = 0) const; + GPUd() void TransformXYZ(int slice, int row, float& x, float& y, float& z, const TPCFastTransform* ref = nullptr, const float* scale = nullptr, int scaleMode = 0) 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; @@ -187,13 +187,13 @@ class TPCFastTransform : public FlatObject 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 TPCFastTransform* ref = nullptr, float scale = 0.f, int scaleMode = 0) const; + GPUd() void InverseTransformYZtoX(int slice, int row, float y, float z, float& x, const TPCFastTransform* ref = nullptr, const float* scale = nullptr, int scaleMode = 0) 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 TPCFastTransform* ref = nullptr, float scale = 0.f, int scaleMode = 0) const; + GPUd() void InverseTransformYZtoNominalYZ(int slice, int row, float y, float z, float& ny, float& nz, const TPCFastTransform* ref = nullptr, const float* scale = nullptr, int scaleMode = 0) 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 TPCFastTransform* ref = nullptr, float scale = 0.f, int scaleMode = 0) const; + GPUd() void InverseTransformXYZtoNominalXYZ(int slice, int row, float x, float y, float z, float& nx, float& ny, float& nz, const TPCFastTransform* ref = nullptr, const float* scale = nullptr, int scaleMode = 0) 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; @@ -262,6 +262,104 @@ class TPCFastTransform : public FlatObject /// maximal possible drift time of the active area GPUd() float getMaxDriftTime(int slice) const; + GPUd() float getScaling(int slice, float z, const float* scale) const + { + const int idxValOffs = (slice < getGeometry().getNumberOfSlicesA()) ? 0 : NLUMIZBINS; + + float scaleVal = 0; + // find closest z bins + const float absZ = GPUCommonMath::Abs(z); + const int offsZIdx = 2 * NLUMIZBINS; + const float zMin = scale[offsZIdx]; + const float zMax = scale[offsZIdx + NLUMIZBINS - 1]; + + int type = 0; + if (absZ < zMin) { + type = 0; + // if z value is smaller than min value return last valid value + scaleVal = scale[idxValOffs]; + } else if (absZ > zMax) { + type = 1; + // if z value is larger than max value return last valid value + scaleVal = scale[idxValOffs + NLUMIZBINS - 1]; + } else { + // if z value is in between, perform linear interpolation + + int idxZClosest = 0; + int idxZClosestN = 1; + float deltaZClosest = GPUCommonMath::Abs(scale[offsZIdx] - z); + + for (int iz = 1; iz < NLUMIZBINS; ++iz) { + const int idxZ = iz + offsZIdx; + const float deltaZScaler = GPUCommonMath::Abs(scale[idxZ] - z); + if (deltaZScaler < deltaZClosest) { + deltaZClosest = deltaZScaler; + idxZClosest = iz; + if (absZ == scale[idxZ]) { + idxZClosestN = iz; + } else if (absZ < scale[idxZ]) { + idxZClosestN = iz - 1; + } else { + idxZClosestN = iz + 1; + } + } + } + if (idxZClosest != idxZClosestN) { + type = 2; + // interpolate linearly now + int idx0 = idxZClosest; + int idx1 = idxZClosestN; + // if left index is larger than right, swap indices + if (idx0 > idx1) { + idx0 = idxZClosestN; + idx1 = idxZClosest; + } + const float y0 = scale[idxValOffs + idx0]; + const float y1 = scale[idxValOffs + idx1]; + const float x0 = scale[offsZIdx + idx0]; + const float x1 = scale[offsZIdx + idx1]; + const float x = absZ; + scaleVal = (y0 * (x1 - x) + y1 * (x - x0)) / (x1 - x0); + + GPUCA_DEBUG_STREAMER_CHECK(if (o2::utils::DebugStreamer::checkStream(o2::utils::StreamFlags::streamFastTransform)) { + std::vector scaler(scale, scale + 3 * NLUMIZBINS); + o2::utils::DebugStreamer::instance()->getStreamer("debug_fasttransform", "UPDATE") << o2::utils::DebugStreamer::instance()->getUniqueTreeName("tree_getScaling_inter").data() + << "slice=" << slice + << "scaleVal=" << scaleVal + << "scaler=" << scaler + << "z=" << z + << "idx0=" << idx0 + << "idx1=" << idx1 + << "y0=" << y0 + << "y1=" << y1 + << "x0=" << x0 + << "x1=" << x1 + << "idxValOffs=" << idxValOffs + << "deltaZClosest=" << deltaZClosest + << "type=" << type + << "\n"; + }) + } else { + type = 3; + // no interpolation required + scaleVal = scale[idxValOffs + idxZClosest]; + } + } + + GPUCA_DEBUG_STREAMER_CHECK(if (o2::utils::DebugStreamer::checkStream(o2::utils::StreamFlags::streamFastTransform)) { + std::vector scaler(scale, scale + 3 * NLUMIZBINS); + o2::utils::DebugStreamer::instance()->getStreamer("debug_fasttransform", "UPDATE") << o2::utils::DebugStreamer::instance()->getUniqueTreeName("tree_getScaling").data() + << "slice=" << slice + << "scaleVal=" << scaleVal + << "scaler=" << scaler + << "z=" << z + << "idxValOffs=" << idxValOffs + << "type=" << type + << "\n"; + }) + return scaleVal; + } + #if !defined(GPUCA_GPUCODE) && !defined(GPUCA_STANDALONE) && !defined(GPUCA_ALIROOT_LIB) int writeToFile(std::string outFName = "", std::string name = ""); @@ -274,7 +372,10 @@ class TPCFastTransform : public FlatObject void setSlowTPCSCCorrection(TFile& inpf); /// \return returns the space charge object which is used for the slow correction - const auto& getCorrectionSlow() const { return *mCorrectionSlow; } + const auto& getCorrectionSlow() const + { + return *mCorrectionSlow; + } #endif // !GPUCA_GPUCODE @@ -329,14 +430,15 @@ class TPCFastTransform : public FlatObject float mLumi; ///< luminosity estimator float mLumiError; ///< error on luminosity float mLumiScaleFactor; ///< user correction factor for lumi (e.g. normalization, efficiency correction etc.) + constexpr static int NLUMIZBINS = 10; /// Correction of (x,u,v) with tricubic interpolator on a regular grid TPCSlowSpaceChargeCorrection* mCorrectionSlow{nullptr}; ///< reference space charge corrections - GPUd() void TransformInternal(int slice, int row, float& u, float& v, float& x, const TPCFastTransform* ref, float scale, int scaleMode) const; + GPUd() void TransformInternal(int slice, int row, float& u, float& v, float& x, const TPCFastTransform* ref, const float* scale, int scaleMode) const; #ifndef GPUCA_ALIROOT_LIB - ClassDefNV(TPCFastTransform, 3); + ClassDefNV(TPCFastTransform, 4); #endif }; @@ -438,11 +540,11 @@ GPUdi() void TPCFastTransform::getTOFcorrection(int slice, int /*row*/, float x, dz = sideC ? dv : -dv; } -GPUdi() void TPCFastTransform::TransformInternal(int slice, int row, float& u, float& v, float& x, const TPCFastTransform* ref, float scale, int scaleMode) const +GPUdi() void TPCFastTransform::TransformInternal(int slice, int row, float& u, float& v, float& x, const TPCFastTransform* ref, const float* scale, int scaleMode) const { if (mApplyCorrection) { float dx = 0.f, du = 0.f, dv = 0.f; - if ((scale >= 0.f) || (scaleMode == 1) || (scaleMode == 2)) { + if ((scale != nullptr) || (scaleMode == 1) || (scaleMode == 2)) { #ifndef GPUCA_GPUCODE if (mCorrectionSlow) { float ly, lz; @@ -464,23 +566,26 @@ GPUdi() void TPCFastTransform::TransformInternal(int slice, int row, float& u, f { mCorrection.getCorrection(slice, row, u, v, dx, du, dv); if (ref) { - if ((scale > 0.f) && (scaleMode == 0)) { // scaling was requested + float scaleVal = getScaling(slice, v, scale); + if ((scaleVal > 0.f) && (scaleMode == 0)) { // scaling was requested float dxRef, duRef, dvRef; ref->mCorrection.getCorrection(slice, row, u, v, dxRef, duRef, dvRef); - dx = (dx - dxRef) * scale + dxRef; - du = (du - duRef) * scale + duRef; - dv = (dv - dvRef) * scale + dvRef; - } else if ((scale != 0.f) && ((scaleMode == 1) || (scaleMode == 2))) { + dx = (dx - dxRef) * scaleVal + dxRef; + du = (du - duRef) * scaleVal + duRef; + dv = (dv - dvRef) * scaleVal + dvRef; + } else if ((scaleVal != 0.f) && ((scaleMode == 1) || (scaleMode == 2))) { float dxRef, duRef, dvRef; ref->mCorrection.getCorrection(slice, row, u, v, dxRef, duRef, dvRef); - dx = dxRef * scale + dx; - du = duRef * scale + du; - dv = dvRef * scale + dv; + dx = dxRef * scaleVal + dx; + du = duRef * scaleVal + du; + dv = dvRef * scaleVal + dv; } } } } GPUCA_DEBUG_STREAMER_CHECK(if (o2::utils::DebugStreamer::checkStream(o2::utils::StreamFlags::streamFastTransform)) { + float scaleVal = getScaling(slice, v, scale); + float ly, lz; getGeometry().convUVtoLocal(slice, u, v, ly, lz); @@ -521,7 +626,7 @@ GPUdi() void TPCFastTransform::TransformInternal(int slice, int row, float& u, f << "u=" << u << "row=" << row << "slice=" << slice - << "scale=" << scale + << "scale=" << scaleVal // original local coordinates << "ly=" << ly << "lz=" << lz @@ -548,7 +653,7 @@ GPUdi() void TPCFastTransform::TransformInternal(int slice, int row, float& u, f } } -GPUdi() void TPCFastTransform::TransformXYZ(int slice, int row, float& x, float& y, float& z, const TPCFastTransform* ref, float scale, int scaleMode) const +GPUdi() void TPCFastTransform::TransformXYZ(int slice, int row, float& x, float& y, float& z, const TPCFastTransform* ref, const float* scale, int scaleMode) const { float u, v; getGeometry().convLocalToUV(slice, y, z, u, v); @@ -559,7 +664,7 @@ GPUdi() void TPCFastTransform::TransformXYZ(int slice, int row, float& x, float& z += dzTOF; } -GPUdi() void TPCFastTransform::Transform(int slice, int row, float pad, float time, float& x, float& y, float& z, float vertexTime, const TPCFastTransform* ref, float scale, int scaleMode) const +GPUdi() void TPCFastTransform::Transform(int slice, int row, float pad, float time, float& x, float& y, float& z, float vertexTime, const TPCFastTransform* ref, const float* scale, int scaleMode) const { /// _______________ The main method: cluster transformation _______________________ /// @@ -740,32 +845,34 @@ GPUdi() float TPCFastTransform::getMaxDriftTime(int slice) const return maxTime; } -GPUdi() void TPCFastTransform::InverseTransformYZtoX(int slice, int row, float y, float z, float& x, const TPCFastTransform* ref, float scale, int scaleMode) const +GPUdi() void TPCFastTransform::InverseTransformYZtoX(int slice, int row, float y, float z, float& x, const TPCFastTransform* ref, const float* scale, int scaleMode) const { /// Transformation y,z -> x float u = 0, v = 0; getGeometry().convLocalToUV(slice, y, z, u, v); - if ((scale >= 0.f) || (scaleMode == 1) || (scaleMode == 2)) { + if ((scale != nullptr) || (scaleMode == 1) || (scaleMode == 2)) { mCorrection.getCorrectionInvCorrectedX(slice, row, u, v, x); if (ref) { // scaling was requested - if (scaleMode == 0 && scale > 0.f) { + float scaleVal = getScaling(slice, v, scale); + if (scaleMode == 0 && scaleVal > 0.f) { float xr; ref->mCorrection.getCorrectionInvCorrectedX(slice, row, u, v, xr); - x = (x - xr) * scale + xr; - } else if ((scale != 0) && ((scaleMode == 1) || (scaleMode == 2))) { + x = (x - xr) * scaleVal + xr; + } else if ((scaleVal != 0) && ((scaleMode == 1) || (scaleMode == 2))) { float xr; ref->mCorrection.getCorrectionInvCorrectedX(slice, row, u, v, xr); - x = (xr - getGeometry().getRowInfo(row).x) * scale + x; // xr=mGeo.getRowInfo(row).x + dx; + x = (xr - getGeometry().getRowInfo(row).x) * scaleVal + x; // xr=mGeo.getRowInfo(row).x + dx; } } } else { x = mCorrection.getGeometry().getRowInfo(row).x; // corrections are disabled } GPUCA_DEBUG_STREAMER_CHECK(if (o2::utils::DebugStreamer::checkStream(o2::utils::StreamFlags::streamFastTransform)) { + float scaleVal = getScaling(slice, v, scale); o2::utils::DebugStreamer::instance()->getStreamer("debug_fasttransform", "UPDATE") << o2::utils::DebugStreamer::instance()->getUniqueTreeName("tree_InverseTransformYZtoX").data() << "slice=" << slice << "row=" << row - << "scale=" << scale + << "scale=" << scaleVal << "y=" << y << "z=" << z << "x=" << x @@ -775,24 +882,25 @@ GPUdi() void TPCFastTransform::InverseTransformYZtoX(int slice, int row, float y }) } -GPUdi() void TPCFastTransform::InverseTransformYZtoNominalYZ(int slice, int row, float y, float z, float& ny, float& nz, const TPCFastTransform* ref, float scale, int scaleMode) const +GPUdi() void TPCFastTransform::InverseTransformYZtoNominalYZ(int slice, int row, float y, float z, float& ny, float& nz, const TPCFastTransform* ref, const float* scale, int scaleMode) const { /// Transformation y,z -> x float u = 0, v = 0, un = 0, vn = 0; getGeometry().convLocalToUV(slice, y, z, u, v); - if ((scale >= 0.f) || (scaleMode == 1) || (scaleMode == 2)) { + if ((scale != nullptr) || (scaleMode == 1) || (scaleMode == 2)) { mCorrection.getCorrectionInvUV(slice, row, u, v, un, vn); if (ref) { // scaling was requested - if (scaleMode == 0 && scale > 0.f) { + float scaleVal = getScaling(slice, v, scale); + if (scaleMode == 0 && scaleVal > 0.f) { float unr = 0, vnr = 0; ref->mCorrection.getCorrectionInvUV(slice, row, u, v, unr, vnr); - un = (un - unr) * scale + unr; - vn = (vn - vnr) * scale + vnr; - } else if ((scale != 0) && ((scaleMode == 1) || (scaleMode == 2))) { + un = (un - unr) * scaleVal + unr; + vn = (vn - vnr) * scaleVal + vnr; + } else if ((scaleVal != 0) && ((scaleMode == 1) || (scaleMode == 2))) { float unr = 0, vnr = 0; ref->mCorrection.getCorrectionInvUV(slice, row, u, v, unr, vnr); - un = (unr - u) * scale + un; // unr = u - duv[0]; - vn = (vnr - v) * scale + vn; + un = (unr - u) * scaleVal + un; // unr = u - duv[0]; + vn = (vnr - v) * scaleVal + vn; } } } else { @@ -802,10 +910,11 @@ GPUdi() void TPCFastTransform::InverseTransformYZtoNominalYZ(int slice, int row, getGeometry().convUVtoLocal(slice, un, vn, ny, nz); GPUCA_DEBUG_STREAMER_CHECK(if (o2::utils::DebugStreamer::checkStream(o2::utils::StreamFlags::streamFastTransform)) { + float scaleVal = getScaling(slice, v, scale); o2::utils::DebugStreamer::instance()->getStreamer("debug_fasttransform", "UPDATE") << o2::utils::DebugStreamer::instance()->getUniqueTreeName("tree_InverseTransformYZtoNominalYZ").data() << "slice=" << slice << "row=" << row - << "scale=" << scale + << "scale=" << scaleVal << "y=" << y << "z=" << z << "ny=" << ny @@ -818,7 +927,7 @@ GPUdi() void TPCFastTransform::InverseTransformYZtoNominalYZ(int slice, int row, }) } -GPUdi() void TPCFastTransform::InverseTransformXYZtoNominalXYZ(int slice, int row, float x, float y, float z, float& nx, float& ny, float& nz, const TPCFastTransform* ref, float scale, int scaleMode) const +GPUdi() void TPCFastTransform::InverseTransformXYZtoNominalXYZ(int slice, int row, float x, float y, float z, float& nx, float& ny, float& nz, const TPCFastTransform* ref, const float* scale, int scaleMode) const { /// Inverse transformation: Transformed X, Y and Z -> X, Y and Z, transformed w/o space charge correction int row2 = row + 1;