Skip to content

Commit

Permalink
Adding z-dependent scalers for 10 z bins
Browse files Browse the repository at this point in the history
  • Loading branch information
@matthias-kleiner
matthias-kleiner committed Jan 3, 2024
1 parent 30b497b commit 62c2e49
Show file tree
Hide file tree
Showing 7 changed files with 299 additions and 106 deletions.
49 changes: 30 additions & 19 deletions Detectors/TPC/calibration/include/TPCCalibration/TPCScaler.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<float> 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<float> mWeightsA; ///< stored weights for A-side
std::vector<float> mWeightsC; ///< stored weights for C-side
std::vector<float> mZPos; ///< z positions

ClassDefNV(TPCScalerWeights, 1);
};
Expand Down Expand Up @@ -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<float> 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); }
Expand All @@ -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<float> 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<float> mScalerA{}; ///< TPC scaler for A-side
std::vector<float> 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<float> mScalerA{}; ///< TPC scaler for A-side
std::vector<float> 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; }
Expand Down
17 changes: 13 additions & 4 deletions Detectors/TPC/calibration/src/CorrectionMapsLoader.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<float>("tpcscaler");
setInstLumi(mInstLumiFactor * tpcScaler);
std::vector<float> tpcScaler = pc.inputs().get<std::vector<float>>("tpcscaler");

if (tpcScaler.size() != getNBinsLumi()) {
const float v[getNBinsLumi()]{mInstLumiFactor * tpcScaler[0], mInstLumiFactor * tpcScaler[1]};
setInstLumi(v);
} else {
setInstLumi(tpcScaler.data());
}
}
}

Expand Down Expand Up @@ -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();
Expand Down
58 changes: 37 additions & 21 deletions Detectors/TPC/calibration/src/TPCScaler.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@
#include <TTree.h>
#include "Framework/Logger.h"
#include "CommonConstants/LHCConstants.h"
#include <numeric>

using namespace o2::tpc;

Expand Down Expand Up @@ -113,38 +114,53 @@ void TPCScaler::setFromTree(TTree& tpcScalerTree)
tpcScalerTree.SetBranchAddress("TPCScaler", nullptr);
}

float TPCScaler::getMeanScaler(double timestamp, o2::tpc::Side side) const
std::vector<float> 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<float>(-1);
LOGP(error, "Empty data provided {}", nValues);
}

// clamp indices to min and max
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<float> 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)
Expand All @@ -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)];
}
41 changes: 31 additions & 10 deletions Detectors/TPC/workflow/src/TPCScalerSpec.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -62,11 +62,29 @@ class TPCScalerSpec : public Task
const auto orbitResetTimeMS = o2::base::GRPGeomHelper::instance().getOrbitResetTimeMS();
const auto firstTFOrbit = pc.services().get<o2::framework::TimingInfo>().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<float> scalerA{mTPCScaler.getMeanScaler(timestamp, o2::tpc::Side::A)};
std::vector<float> scalerC{mTPCScaler.getMeanScaler(timestamp, o2::tpc::Side::C)};
const std::vector<float> 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
Expand All @@ -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);
}
}
Expand All @@ -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)
Expand Down
6 changes: 4 additions & 2 deletions GPU/TPCFastTransformation/CorrectionMapsHelper.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
}
Expand Down Expand Up @@ -78,5 +80,5 @@ void CorrectionMapsHelper::setCorrMapRef(std::unique_ptr<TPCFastTransform>&& 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());
}
49 changes: 37 additions & 12 deletions GPU/TPCFastTransformation/CorrectionMapsHelper.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
}
}
Expand Down Expand Up @@ -97,25 +103,39 @@ 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) {
reportScaling();
}
}

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; }
Expand Down Expand Up @@ -146,17 +166,22 @@ 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,
LumiBit = 0x4 };
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
Expand Down
Loading

0 comments on commit 62c2e49

Please sign in to comment.