From ab1336c3f6a2e7ba4a5f69483be83afe04fb3700 Mon Sep 17 00:00:00 2001 From: Matthias Kleiner Date: Fri, 29 Dec 2023 09:12:11 +0100 Subject: [PATCH] Adding TPC scaler weights for 1D-distortion fluctuation correction --- Detectors/TPC/base/include/TPCBase/CDBTypes.h | 2 + .../include/TPCCalibration/TPCScaler.h | 66 +++++++++++-- .../calibration/src/TPCCalibrationLinkDef.h | 1 + Detectors/TPC/calibration/src/TPCScaler.cxx | 98 ++++++++++++++++++- .../include/TPCWorkflow/TPCScalerSpec.h | 2 +- Detectors/TPC/workflow/src/TPCScalerSpec.cxx | 30 +++++- Detectors/TPC/workflow/src/tpc-scaler.cxx | 7 +- 7 files changed, 187 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..03a55a4b732b4 100644 --- a/Detectors/TPC/calibration/include/TPCCalibration/TPCScaler.h +++ b/Detectors/TPC/calibration/include/TPCCalibration/TPCScaler.h @@ -27,6 +27,16 @@ 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 mSamplingTimeMS = -1; ///< sampling of the stored weights + float mFirstTimeStampMS = 0; ///< first timestamp + std::vector mWeights; ///< stored weights + + ClassDefNV(TPCScalerWeights, 1); +}; + class TPCScaler { public: @@ -62,6 +72,19 @@ class TPCScaler /// \param file output file void dumpToFile(const char* file, const char* name); + /// dump this object to a file for a given time range + /// \param file output file + /// \param startTime start time stamp in ms of the dumped object + /// \param endTime end time stamp in ms of the dumped object (if endTimeMS=-1 the end of the stored data is used) + void dumpToFile(const char* file, const char* name, double startTimeMS, double endTimeMS); + + /// dump this object to several file each containing the scalers for n minutes + /// \param file output file + /// \param minutesPerObject minutes each output file is covering + /// \param marginMS additional margin to the beginning of each object (should be >= getIonDriftTimeMS()) + /// \param marginCCDBMinutes margin for CCDB timestamp + void dumpToFileSlices(const char* file, const char* name, int minutesPerObject, float marginMS = 500, float marginCCDBMinutes = 1); + /// load parameters from input file (which were written using the writeToFile method) /// \param inpf input file void loadFromFile(const char* inpf, const char* name); @@ -87,6 +110,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 +123,38 @@ 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); } + + /// clamp scaler values + /// \param minThreshold minimum accepted scaler value + /// \param maxThreshold maximum accepted scaler value + 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; } + + /// \return returns stored weights for TPC scalers + const auto& getScalerWeights(o2::tpc::Side side) const { return (side == o2::tpc::Side::A) ? mScalerWeightsA : mScalerWeightsC; } + + void useWeights(bool useWeights) { mUseWeights = useWeights; } + 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 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 + + // 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..d705232bfdd91 100644 --- a/Detectors/TPC/calibration/src/TPCScaler.cxx +++ b/Detectors/TPC/calibration/src/TPCScaler.cxx @@ -18,18 +18,81 @@ #include #include #include "Framework/Logger.h" +#include "CommonConstants/LHCConstants.h" 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); } +void TPCScaler::dumpToFile(const char* file, const char* name, double startTimeMS, double endTimeMS) +{ + TPCScaler scaler; + scaler.mIonDriftTimeMS = mIonDriftTimeMS; + scaler.mRun = mRun; + scaler.mFirstTFOrbit = (startTimeMS <= 0) ? mFirstTFOrbit : (mFirstTFOrbit + (startTimeMS - mTimeStampMS) / (o2::constants::lhc::LHCOrbitMUS * 0.001)); + scaler.mTimeStampMS = (startTimeMS <= 0) ? mTimeStampMS : startTimeMS; + scaler.mIntegrationTimeMS = mIntegrationTimeMS; + + const int dataIdx = (startTimeMS <= 0) ? 0 : getDataIdx(startTimeMS); + LOGP(info, "dataIdx: {} max: {}", dataIdx, getNValues(Side::A)); + const int idxDataStart = (startTimeMS <= 0) ? 0 : dataIdx; + const int idxDataEndA = ((endTimeMS < 0) || (dataIdx > getNValues(Side::A))) ? getNValues(Side::A) : dataIdx; + const int idxDataEndC = ((endTimeMS < 0) || (dataIdx > getNValues(Side::C))) ? getNValues(Side::C) : dataIdx; + scaler.mScalerA = std::vector(mScalerA.begin() + idxDataStart, mScalerA.begin() + idxDataEndA); + scaler.mScalerC = std::vector(mScalerC.begin() + idxDataStart, mScalerC.begin() + idxDataEndC); + scaler.dumpToFile(file, name); +} + +void TPCScaler::dumpToFileSlices(const char* file, const char* name, int minutesPerObject, float marginMS, float marginCCDBMinutes) +{ + if (getNValues(o2::tpc::Side::A) != getNValues(o2::tpc::Side::C)) { + LOGP(error, "Number of points stored for A-side and C-side is different"); + return; + } + const long marginCCDBMS = marginCCDBMinutes * 60 * 1000; + const float msPerObjectTmp = minutesPerObject * 60 * 1000; + const int valuesPerSlice = static_cast(msPerObjectTmp / mIntegrationTimeMS); + const int marginPerSlice = static_cast(marginMS / mIntegrationTimeMS); + int nSlices = getNValues(Side::A) / valuesPerSlice; // number of output objects + + LOGP(info, "Producing {} objects with a CCDB margin of {} ms and {} margin per slice", nSlices, marginCCDBMS, marginPerSlice); + if (nSlices == 0) { + nSlices = 1; + } + + for (int i = 0; i < nSlices; ++i) { + const int idxDataStart = (i == 0) ? 0 : (i * valuesPerSlice - marginPerSlice); + int idxDataEnd = (i == nSlices - 1) ? getNValues(Side::A) : ((i + 1) * valuesPerSlice + marginPerSlice); + if (idxDataEnd > getNValues(Side::A)) { + idxDataEnd = getNValues(Side::A); + } + + TPCScaler scaler; + scaler.mIonDriftTimeMS = mIonDriftTimeMS; + scaler.mRun = mRun; + scaler.mTimeStampMS = (i == 0) ? mTimeStampMS : (mTimeStampMS + idxDataStart * static_cast(mIntegrationTimeMS)); + scaler.mFirstTFOrbit = mFirstTFOrbit + (scaler.mTimeStampMS - mTimeStampMS) / (o2::constants::lhc::LHCOrbitMUS * 0.001); + scaler.mIntegrationTimeMS = mIntegrationTimeMS; + scaler.mScalerA = std::vector(mScalerA.begin() + idxDataStart, mScalerA.begin() + idxDataEnd); + scaler.mScalerC = std::vector(mScalerC.begin() + idxDataStart, mScalerC.begin() + idxDataEnd); + + const long timePerSliceMS = valuesPerSlice * mIntegrationTimeMS; + const long tsCCDBStart = mTimeStampMS + i * timePerSliceMS; + const long tsCCDBStartMargin = (i == 0) ? (tsCCDBStart - marginCCDBMS) : tsCCDBStart; + const long tsCCDBEnd = (i == nSlices - 1) ? (getEndTimeStampMS(o2::tpc::Side::A) + marginCCDBMS) : (tsCCDBStart + timePerSliceMS); + const std::string fileOut = fmt::format("{}_{}_{}_{}.root", file, i, tsCCDBStartMargin, tsCCDBEnd - 1); + scaler.dumpToFile(fileOut.data(), name); + } +} + void TPCScaler::loadFromFile(const char* inpf, const char* name) { TFile out(inpf, "READ"); @@ -53,7 +116,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 +130,35 @@ 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(side).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(side).getWeight(relTSMS); + } + sum += getScalers(i, side) * weight; + sumW += weight; } - return (sum / nVals); + if (sumW != 0) { + return (sum / sumW); + } + return 0; +} + +void TPCScaler::clampScalers(float minThreshold, float maxThreshold, Side side) +{ + auto& scaler = (side == o2::tpc::Side::A) ? mScalerA : mScalerC; + 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 +{ + const int idx = (deltaTime - mFirstTimeStampMS) / mSamplingTimeMS + 0.5; + if ((idx < 0) || (idx > mWeights.size() - 1)) { + // LOGP(info, "Index out of range for deltaTime: {} mFirstTimeStampMS: {} mSamplingTimeMS: {}", deltaTime, mFirstTimeStampMS, mSamplingTimeMS); + return 1; + } + return mWeights[idx]; +} \ No newline at end of file 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..4f23623532f3f 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,37 @@ 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, Side::A); + mTPCScaler.setScalerWeights(mScalerWeights, Side::C); + mTPCScaler.useWeights(true); + } + } + if (matcher == ConcreteDataMatcher(o2::header::gDataOriginTPC, "TPCSCALERCCDBW", 0)) { + LOGP(info, "Updating TPC scaler weights"); + mScalerWeights = *(TPCScalerWeights*)obj; + mTPCScaler.setScalerWeights(mScalerWeights, Side::A); + mTPCScaler.setScalerWeights(mScalerWeights, Side::C); + 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 +126,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; }