From deba975d43f3ee2478c083c1fa683f0d842a44b9 Mon Sep 17 00:00:00 2001 From: Matthias Kleiner Date: Fri, 29 Dec 2023 09:12:11 +0100 Subject: [PATCH] TPC: Adding scaler weights for 1D-distortion fluctuation correction --- Detectors/TPC/base/include/TPCBase/CDBTypes.h | 2 + .../include/TPCCalibration/TPCScaler.h | 52 +++++++++++++++---- .../calibration/src/TPCCalibrationLinkDef.h | 1 + Detectors/TPC/calibration/src/TPCScaler.cxx | 45 ++++++++++++++-- .../include/TPCWorkflow/TPCScalerSpec.h | 2 +- Detectors/TPC/workflow/src/TPCScalerSpec.cxx | 28 ++++++++-- Detectors/TPC/workflow/src/tpc-scaler.cxx | 7 ++- 7 files changed, 118 insertions(+), 19 deletions(-) diff --git a/Detectors/TPC/base/include/TPCBase/CDBTypes.h b/Detectors/TPC/base/include/TPCBase/CDBTypes.h index 78ffacbd552fa..48a79b36d96b8 100644 --- a/Detectors/TPC/base/include/TPCBase/CDBTypes.h +++ b/Detectors/TPC/base/include/TPCBase/CDBTypes.h @@ -78,6 +78,7 @@ enum class CDBType { /// CalTimeSeries, ///< integrated DCAs for longer time interval CalScaler, ///< Scaler from IDCs or combined estimator + CalScalerWeights, ///< Weights for scalers /// CorrMapParam, ///< parameters for CorrectionMapsLoader configuration /// @@ -142,6 +143,7 @@ const std::unordered_map CDBTypeMap{ // time series {CDBType::CalTimeSeries, "TPC/Calib/TimeSeries"}, {CDBType::CalScaler, "TPC/Calib/Scaler"}, + {CDBType::CalScalerWeights, "TPC/Calib/ScalerWeights"}, // correction maps loader params {CDBType::CorrMapParam, "TPC/Calib/CorrMapParam"}, // distortion maps diff --git a/Detectors/TPC/calibration/include/TPCCalibration/TPCScaler.h b/Detectors/TPC/calibration/include/TPCCalibration/TPCScaler.h index 131fb626c7654..532c263244394 100644 --- a/Detectors/TPC/calibration/include/TPCCalibration/TPCScaler.h +++ b/Detectors/TPC/calibration/include/TPCCalibration/TPCScaler.h @@ -27,6 +27,17 @@ namespace o2::tpc Class for storing the scalers which are used to calculate an estimate for the mean space-charge density for the last ion drift time */ +struct TPCScalerWeights { + float getWeight(float deltaTime) const; + bool isValid() const { return (mSamplingTimeMS > 0 && !mWeights.empty()); } + float getDurationMS() const { return mSamplingTimeMS * mWeights.size(); } + float mSamplingTimeMS = -1; ///< sampling of the stored weights + float mFirstTimeStampMS = 0; ///< first timestamp + std::vector mWeights; ///< stored weights + + ClassDefNV(TPCScalerWeights, 1); +}; + class TPCScaler { public: @@ -87,6 +98,9 @@ class TPCScaler /// \return return first time in ms double getStartTimeStampMS() const { return mTimeStampMS; } + /// \return return end time in ms + double getEndTimeStampMS(o2::tpc::Side side) const { return mTimeStampMS + getDurationMS(side); } + /// \return returns integration time for each scaler value float getIntegrationTimeMS() const { return mIntegrationTimeMS; } @@ -97,16 +111,36 @@ class TPCScaler /// \param timestamp timestamp for which the last values are used to calculate the mean float getMeanScaler(double timestamp, 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); } + + /// setting the weights for the scalers + void setScalerWeights(const TPCScalerWeights& weights) { mScalerWeights = weights; } + + /// \return returns stored weights for TPC scalers + const auto& getScalerWeights() const { return mScalerWeights; } + + /// enable usage of weights + void useWeights(bool useWeights) { mUseWeights = useWeights; } + + /// return if weights are used + bool weightsUsed() const { return mUseWeights; } + 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 - - ClassDefNV(TPCScaler, 1); + 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; } + + ClassDefNV(TPCScaler, 2); }; } // namespace o2::tpc diff --git a/Detectors/TPC/calibration/src/TPCCalibrationLinkDef.h b/Detectors/TPC/calibration/src/TPCCalibrationLinkDef.h index 238482b6f8687..9b0fd7ce71e0e 100644 --- a/Detectors/TPC/calibration/src/TPCCalibrationLinkDef.h +++ b/Detectors/TPC/calibration/src/TPCCalibrationLinkDef.h @@ -113,4 +113,5 @@ #pragma link C++ class o2::tpc::CalculatedEdx + ; #pragma link C++ class o2::tpc::TPCScaler + ; +#pragma link C++ struct o2::tpc::TPCScalerWeights + ; #endif diff --git a/Detectors/TPC/calibration/src/TPCScaler.cxx b/Detectors/TPC/calibration/src/TPCScaler.cxx index cb145dd112b07..f10555892bb33 100644 --- a/Detectors/TPC/calibration/src/TPCScaler.cxx +++ b/Detectors/TPC/calibration/src/TPCScaler.cxx @@ -24,7 +24,8 @@ using namespace o2::tpc; void TPCScaler::dumpToFile(const char* file, const char* name) { TFile out(file, "RECREATE"); - TTree tree("TPCScaler", "TPCScaler"); + TTree tree(name, name); + tree.SetAutoSave(0); tree.Branch("TPCScaler", this); tree.Fill(); out.WriteObject(&tree, name); @@ -53,7 +54,7 @@ void TPCScaler::setFromTree(TTree& tpcScalerTree) float TPCScaler::getMeanScaler(double timestamp, o2::tpc::Side side) const { // index to data buffer - const int idxData = (timestamp - mTimeStampMS) / mIntegrationTimeMS + 0.5; + const int idxData = getDataIdx(timestamp); const int nVals = getNValuesIonDriftTime(); const int nValues = getNValues(side); if ((nVals == 0) || (nVals > nValues)) { @@ -67,8 +68,44 @@ float TPCScaler::getMeanScaler(double timestamp, o2::tpc::Side side) const // sump up values from last ion drift time float sum = 0; + float sumW = 0; + const bool useWeights = mUseWeights && getScalerWeights().isValid(); for (int i = firstIdx; i < lastIdx; ++i) { - sum += getScalers(i, side); + float weight = 1; + if (useWeights) { + const double relTSMS = mTimeStampMS + i * mIntegrationTimeMS - timestamp; + weight = getScalerWeights().getWeight(relTSMS); + } + sum += getScalers(i, side) * weight; + sumW += weight; + } + if (sumW != 0) { + return (sum / sumW); + } + return 0; +} + +float TPCScalerWeights::getWeight(float deltaTime) const +{ + const float idxF = (deltaTime - mFirstTimeStampMS) / mSamplingTimeMS; + const int idx = idxF; + if ((idx < 0) || (idx > mWeights.size() - 1)) { + LOGP(error, "Index out of range for deltaTime: {} mFirstTimeStampMS: {} mSamplingTimeMS: {}", deltaTime, mFirstTimeStampMS, mSamplingTimeMS); + // set weight 1 to in case it is out of range. This can only happen if the TPC scaler is not valid for given time + return 1; + } + + if ((idxF == idx) || (idx == mWeights.size() - 1)) { + // no interpolation required + return mWeights[idx]; + } else { + // interpolate scaler + const float y0 = mWeights[idx]; + const float y1 = mWeights[idx + 1]; + const float x0 = idx; + const float x1 = idx + 1; + const float x = idxF; + const float y = ((y0 * (x1 - x)) + y1 * (x - x0)) / (x1 - x0); + return y; } - return (sum / nVals); } diff --git a/Detectors/TPC/workflow/include/TPCWorkflow/TPCScalerSpec.h b/Detectors/TPC/workflow/include/TPCWorkflow/TPCScalerSpec.h index 2bd3ca89bfb6a..3cea2e0175772 100644 --- a/Detectors/TPC/workflow/include/TPCWorkflow/TPCScalerSpec.h +++ b/Detectors/TPC/workflow/include/TPCWorkflow/TPCScalerSpec.h @@ -19,7 +19,7 @@ namespace o2 namespace tpc { -o2::framework::DataProcessorSpec getTPCScalerSpec(); +o2::framework::DataProcessorSpec getTPCScalerSpec(bool enableWeights); } // end namespace tpc } // end namespace o2 diff --git a/Detectors/TPC/workflow/src/TPCScalerSpec.cxx b/Detectors/TPC/workflow/src/TPCScalerSpec.cxx index 5d33b745a0ae7..198b42f5be565 100644 --- a/Detectors/TPC/workflow/src/TPCScalerSpec.cxx +++ b/Detectors/TPC/workflow/src/TPCScalerSpec.cxx @@ -34,7 +34,7 @@ namespace tpc class TPCScalerSpec : public Task { public: - TPCScalerSpec(std::shared_ptr req) : mCCDBRequest(req){}; + TPCScalerSpec(std::shared_ptr req, bool enableWeights) : mCCDBRequest(req), mEnableWeights(enableWeights){}; void init(framework::InitContext& ic) final { @@ -49,6 +49,12 @@ class TPCScalerSpec : public Task pc.inputs().get("tpcscaler"); } + if (mEnableWeights) { + if (pc.inputs().isValid("tpcscalerw")) { + pc.inputs().get("tpcscalerw"); + } + } + if (pc.services().get().runNumber != mTPCScaler.getRun()) { LOGP(error, "Run number {} of processed data and run number {} of loaded TPC scaler doesnt match!", pc.services().get().runNumber, mTPCScaler.getRun()); } @@ -73,19 +79,35 @@ class TPCScalerSpec : public Task LOGP(info, "Setting ion drift time to: {}", mIonDriftTimeMS); mTPCScaler.setIonDriftTimeMS(mIonDriftTimeMS); } + if (mScalerWeights.isValid()) { + LOGP(info, "Setting TPC scaler weights"); + 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); + mTPCScaler.useWeights(true); } } private: std::shared_ptr mCCDBRequest; ///< info for CCDB request + const bool mEnableWeights{false}; ///< use weights for TPC scalers + TPCScalerWeights mScalerWeights{}; ///< scaler weights float mIonDriftTimeMS{-1}; ///< ion drift time TPCScaler mTPCScaler; ///< tpc scaler }; -o2::framework::DataProcessorSpec getTPCScalerSpec() +o2::framework::DataProcessorSpec getTPCScalerSpec(bool enableWeights) { std::vector inputs; inputs.emplace_back("tpcscaler", o2::header::gDataOriginTPC, "TPCSCALERCCDB", 0, Lifetime::Condition, ccdbParamSpec(o2::tpc::CDBTypeMap.at(o2::tpc::CDBType::CalScaler), {}, 1)); // time-dependent + if (enableWeights) { + inputs.emplace_back("tpcscalerw", o2::header::gDataOriginTPC, "TPCSCALERWCCDB", 0, Lifetime::Condition, ccdbParamSpec(o2::tpc::CDBTypeMap.at(o2::tpc::CDBType::CalScalerWeights), {}, 0)); // non time-dependent + } auto ccdbRequest = std::make_shared(true, // orbitResetTime false, // GRPECS=true for nHBF per TF @@ -102,7 +124,7 @@ o2::framework::DataProcessorSpec getTPCScalerSpec() "tpc-scaler", inputs, outputs, - AlgorithmSpec{adaptFromTask(ccdbRequest)}, + AlgorithmSpec{adaptFromTask(ccdbRequest, enableWeights)}, Options{ {"ion-drift-time", VariantType::Float, -1.f, {"Overwrite ion drift time if a value >0 is provided"}}}}; } diff --git a/Detectors/TPC/workflow/src/tpc-scaler.cxx b/Detectors/TPC/workflow/src/tpc-scaler.cxx index e9bc8c4a1f035..54034d9e04cbb 100644 --- a/Detectors/TPC/workflow/src/tpc-scaler.cxx +++ b/Detectors/TPC/workflow/src/tpc-scaler.cxx @@ -21,7 +21,9 @@ using namespace o2::framework; void customize(std::vector& workflowOptions) { // option allowing to set parameters - std::vector options{ConfigParamSpec{"configKeyValues", VariantType::String, "", {"Semicolon separated key=value strings"}}}; + std::vector options{ + ConfigParamSpec{"configKeyValues", VariantType::String, "", {"Semicolon separated key=value strings"}}, + {"enableWeights", VariantType::Bool, false, {"Enable weights for TPC scalers"}}}; std::swap(workflowOptions, options); } @@ -31,6 +33,7 @@ WorkflowSpec defineDataProcessing(ConfigContext const& config) { WorkflowSpec workflow; o2::conf::ConfigurableParam::updateFromString(config.options().get("configKeyValues")); - workflow.emplace_back(o2::tpc::getTPCScalerSpec()); + const bool enableWeights = config.options().get("enableWeights"); + workflow.emplace_back(o2::tpc::getTPCScalerSpec(enableWeights)); return workflow; }