diff --git a/Common/MathUtils/include/MathUtils/fit.h b/Common/MathUtils/include/MathUtils/fit.h index 00c39486a4ba0..cd5cb415070d3 100644 --- a/Common/MathUtils/include/MathUtils/fit.h +++ b/Common/MathUtils/include/MathUtils/fit.h @@ -20,6 +20,7 @@ #include #include #include +#include #include "Rtypes.h" #include "TLinearFitter.h" @@ -69,9 +70,9 @@ TFitResultPtr fit(const size_t nBins, const T* arr, const T xMin, const T xMax, // create an empty TFitResult std::shared_ptr tfr(new TFitResult()); // create the fitter from an empty fit result - //std::shared_ptr fitter(new ROOT::Fit::Fitter(std::static_pointer_cast(tfr) ) ); + // std::shared_ptr fitter(new ROOT::Fit::Fitter(std::static_pointer_cast(tfr) ) ); ROOT::Fit::Fitter fitter(tfr); - //ROOT::Fit::FitConfig & fitConfig = fitter->Config(); + // ROOT::Fit::FitConfig & fitConfig = fitter->Config(); const double binWidth = double(xMax - xMin) / double(nBins); @@ -225,8 +226,8 @@ bool medmadGaus(size_t nBins, const T* arr, const T xMin, const T xMax, std::arr /// -1: only one point has been used for the calculation - center of gravity was uesed for calculation /// -4: invalid result!! /// -//template -//Double_t fitGaus(const size_t nBins, const T *arr, const T xMin, const T xMax, std::vector& param); +// template +// Double_t fitGaus(const size_t nBins, const T *arr, const T xMin, const T xMax, std::vector& param); template Double_t fitGaus(const size_t nBins, const T* arr, const T xMin, const T xMax, std::vector& param) { @@ -301,7 +302,7 @@ Double_t fitGaus(const size_t nBins, const T* arr, const T xMin, const T xMax, s Double_t chi2 = 0; if (npoints >= 3) { if (npoints == 3) { - //analytic calculation of the parameters for three points + // analytic calculation of the parameters for three points A.Invert(); TMatrixD res(1, 3); res.Mult(A, b); @@ -334,7 +335,7 @@ Double_t fitGaus(const size_t nBins, const T* arr, const T xMin, const T xMax, s } if (npoints == 2) { - //use center of gravity for 2 points + // use center of gravity for 2 points meanCOG /= sumCOG; rms2COG /= sumCOG; param[0] = max; @@ -524,7 +525,7 @@ R median(std::vector v) auto n = v.size() / 2; nth_element(v.begin(), v.begin() + n, v.end()); auto med = R{v[n]}; - if (!(v.size() & 1)) { //If the set size is even + if (!(v.size() & 1)) { // If the set size is even auto max_it = max_element(v.begin(), v.begin() + n); med = R{(*max_it + med) / 2.0}; } @@ -788,6 +789,169 @@ T MAD2Sigma(int np, T* y) return median * 1.4826; // convert to Gaussian sigma } +/// \return returns the index of the closest timestamps to the left and right of the given timestamp +/// \param timestamps vector of timestamps +/// \param timestamp the timestamp to find the closest timestamps for +template +std::optional> findClosestIndices(const std::vector& timestamps, DataTime timestamp) +{ + if (timestamps.empty()) { + LOGP(warning, "Timestamp vector is empty!"); + return std::nullopt; + } + + if (timestamp <= timestamps.front()) { + return std::pair{0, 0}; + } else if (timestamp >= timestamps.back()) { + return std::pair{timestamps.size() - 1, timestamps.size() - 1}; + } + + const auto it = std::lower_bound(timestamps.begin(), timestamps.end(), timestamp); + const size_t idx = std::distance(timestamps.begin(), it); + const auto prevTimestamp = timestamps[idx - 1]; + const auto nextTimestamp = timestamps[idx]; + return std::pair{(idx - 1), idx}; +} + +struct RollingStats { + RollingStats() = default; + RollingStats(const int nValues) + { + median.resize(nValues); + std.resize(nValues); + nPoints.resize(nValues); + closestDistanceL.resize(nValues); + closestDistanceR.resize(nValues); + } + + std::vector median; ///< median of rolling data + std::vector std; ///< std of rolling data + std::vector nPoints; ///< number of points used for the calculation + std::vector closestDistanceL; ///< distance of closest point to the left + std::vector closestDistanceR; ///< distance of closest point to the right + + ClassDefNV(RollingStats, 1); +}; + +/// \brief calculates the rolling statistics of the input data +/// \return returns the rolling statistics +/// \param timeData times of the input data (assumed to be sorted) +/// \param data values of the input data +/// \param times times for which to calculate the rolling statistics +/// \param deltaMax time range for which the rolling statistics is calculated +/// \param mNthreads number of threads to use for the calculation +/// \param minPoints minimum number of points to use for the calculation of the statistics - otherwise use nearest nClosestPoints points weighted with distance +/// \param nClosestPoints number of closest points in case of number of points in given range is smaller than minPoints +template +RollingStats getRollingStatistics(const DataTimeType& timeData, const DataType& data, const DataTime& times, const double deltaMax, const int mNthreads, const size_t minPoints = 4, const size_t nClosestPoints = 4) +{ + // output statistics + const size_t vecSize = times.size(); + RollingStats stats(vecSize); + + if (!std::is_sorted(timeData.begin(), timeData.end())) { + LOGP(error, "Input data is NOT sorted!"); + return stats; + } + + if (timeData.empty()) { + LOGP(error, "Input data is empty!"); + return stats; + } + + const size_t dataSize = data.size(); + const size_t timeDataSize = timeData.size(); + if (timeDataSize != dataSize) { + LOGP(error, "Input data has different sizes {}!={}", timeDataSize, dataSize); + return stats; + } + + auto myThread = [&](int iThread) { + // data in given time window for median calculation + DataType window; + for (size_t i = iThread; i < vecSize; i += mNthreads) { + const double timeI = times[i]; + + // lower index + const double timeStampLower = timeI - deltaMax; + const auto lower = std::lower_bound(timeData.begin(), timeData.end(), timeStampLower); + size_t idxStart = std::distance(timeData.begin(), lower); + + // upper index + const double timeStampUpper = timeI + deltaMax; + const auto upper = std::lower_bound(timeData.begin(), timeData.end(), timeStampUpper); + size_t idxEnd = std::distance(timeData.begin(), upper); + + // closest data point + if (auto idxClosest = findClosestIndices(timeData, timeI)) { + auto [idxLeft, idxRight] = *idxClosest; + const auto closestL = std::abs(timeData[idxLeft] - timeI); + const auto closestR = std::abs(timeData[idxRight] - timeI); + stats.closestDistanceL[i] = closestL; + stats.closestDistanceR[i] = closestR; + + // if no points are in the range use the n closest points - n from the left and n from the right + const size_t reqSize = idxEnd - idxStart; + if (reqSize < minPoints) { + // calculate weighted average + idxStart = (idxRight > nClosestPoints) ? (idxRight - nClosestPoints) : 0; + idxEnd = std::min(data.size(), idxRight + nClosestPoints); + constexpr float epsilon = 1e-6f; + double weightedSum = 0.0; + double weightTotal = 0.0; + for (size_t j = idxStart; j < idxEnd; ++j) { + const double dist = std::abs(timeI - timeData[j]); + const double weight = 1.0 / (dist + epsilon); + weightedSum += weight * data[j]; + weightTotal += weight; + } + stats.median[i] = (weightTotal > 0.) ? (weightedSum / weightTotal) : 0.0f; + } else { + // calculate statistics + stats.nPoints[i] = reqSize; + + if (idxStart >= data.size()) { + stats.median[i] = data.back(); + continue; + } + + if (reqSize <= 1) { + stats.median[i] = data[idxStart]; + continue; + } + + // calculate median + window.clear(); + if (reqSize > window.capacity()) { + window.reserve(static_cast(reqSize * 1.5)); + } + window.insert(window.end(), data.begin() + idxStart, data.begin() + idxEnd); + const size_t middle = window.size() / 2; + std::nth_element(window.begin(), window.begin() + middle, window.end()); + stats.median[i] = (window.size() % 2 == 1) ? window[middle] : ((window[middle - 1] + window[middle]) / 2.0); + + // calculate the stdev + const float mean = std::accumulate(window.begin(), window.end(), 0.0f) / window.size(); + std::transform(window.begin(), window.end(), window.begin(), [mean](const float val) { return val - mean; }); + const float sqsum = std::inner_product(window.begin(), window.end(), window.begin(), 0.0f); + const float stdev = std::sqrt(sqsum / window.size()); + stats.std[i] = stdev; + } + } + } + }; + + std::vector threads(mNthreads); + for (int i = 0; i < mNthreads; i++) { + threads[i] = std::thread(myThread, i); + } + + for (auto& th : threads) { + th.join(); + } + return stats; +} + } // namespace math_utils } // namespace o2 #endif diff --git a/Common/MathUtils/src/MathUtilsLinkDef.h b/Common/MathUtils/src/MathUtilsLinkDef.h index 6067dd540110c..0b070e537afcd 100644 --- a/Common/MathUtils/src/MathUtilsLinkDef.h +++ b/Common/MathUtils/src/MathUtilsLinkDef.h @@ -46,4 +46,6 @@ #pragma link C++ class o2::math_utils::Legendre1DPolynominal + ; #pragma link C++ class o2::math_utils::Legendre2DPolynominal + ; +#pragma link C++ class o2::math_utils::RollingStats + ; + #endif diff --git a/DataFormats/Detectors/TPC/CMakeLists.txt b/DataFormats/Detectors/TPC/CMakeLists.txt index 81b1d5efad59a..2cc69e16001a6 100644 --- a/DataFormats/Detectors/TPC/CMakeLists.txt +++ b/DataFormats/Detectors/TPC/CMakeLists.txt @@ -35,7 +35,8 @@ o2_add_library( O2::CommonDataFormat O2::Headers O2::DataSampling - O2::Algorithm) + O2::Algorithm + ROOT::Minuit) o2_target_root_dictionary( DataFormatsTPC diff --git a/DataFormats/Detectors/TPC/include/DataFormatsTPC/DCS.h b/DataFormats/Detectors/TPC/include/DataFormatsTPC/DCS.h index 2f9f17f164872..3608fdc390203 100644 --- a/DataFormats/Detectors/TPC/include/DataFormatsTPC/DCS.h +++ b/DataFormats/Detectors/TPC/include/DataFormatsTPC/DCS.h @@ -30,9 +30,13 @@ #include "Framework/Logger.h" #include "DataFormatsTPC/Defs.h" +#include "MathUtils/fit.h" using namespace o2::tpc; +class TLinearFitter; +class TTree; + namespace o2::tpc::dcs { @@ -66,6 +70,19 @@ struct DataPointVector { uint32_t sensorNumber{}; std::vector data; + /// \brief convert data points to a vector of pairs: pair.first -> data and pair.second -> time + auto getPairOfVector() const + { + std::pair, std::vector> pairs; + pairs.first.reserve(data.size()); + pairs.second.reserve(data.size()); + for (const auto& dp : data) { + pairs.first.emplace_back(dp.value); + pairs.second.emplace_back(dp.time); + } + return pairs; + } + void fill(const TimeStampType time, const T& value) { data.emplace_back(DPType{time, value}); } void fill(const DPType& dataPoint) { data.emplace_back(dataPoint); } @@ -169,6 +186,45 @@ const T getAverageValueForTime(const std::vector>& dpVec return (nPoints > 0) ? ret / static_cast(nPoints) : T{}; } +template +dcs::TimeStampType getMinTime(const std::vector>& data, const bool roundToInterval, dcs::TimeStampType fitInterval) +{ + constexpr auto max = std::numeric_limits::max(); + dcs::TimeStampType firstTime = std::numeric_limits::max(); + for (const auto& sensor : data) { + const auto time = sensor.data.size() ? sensor.data.front().time : max; + firstTime = std::min(firstTime, time); + } + + // mFitInterval is is seconds. Round to full amount. + // if e.g. mFitInterval = 5min, then round 10:07:20.510 to 10:05:00.000 + if (roundToInterval) { + firstTime -= (firstTime % fitInterval); + } + + return firstTime; +} + +template +dcs::TimeStampType getMaxTime(const std::vector>& data) +{ + constexpr auto min = 0; + dcs::TimeStampType lastTime = 0; + for (const auto& sensor : data) { + const auto time = sensor.data.size() ? sensor.data.back().time : 0; + lastTime = std::max(lastTime, time); + } + + // mFitInterval is is seconds. Round to full amount. + // if e.g. mFitInterval = 5min, then round 10:07:20.510 to 10:05:00.000 + // TODO: fix this + // if (mRoundToInterval) { + // lastTime -= (lastTime % mFitInterval); + //} + + return lastTime; +} + using RawDPsF = DataPointVector; // using RawDPsI = DataPointVector; @@ -210,6 +266,12 @@ struct Temperature { static constexpr auto& getSensorPosition(const size_t sensor) { return SensorPosition[sensor]; } + /// \brief make fit of the mean temperature and gradients in time intervals + /// \param Side TPC side for which to make the fit + /// \param fitInterval time interval for the fits + /// \param roundToInterval round min time + void fitTemperature(Side side, dcs::TimeStampType fitInterval = 5 * 60 * 1000, const bool roundToInterval = false); + struct Stats { DataType mean{}; ///< average temperature in K DataType gradX{}; ///< horizontal temperature gradient in K/cm @@ -262,6 +324,9 @@ struct Temperature { doAppend(raw, other.raw); } + private: + bool makeFit(TLinearFitter& fitter, const int nDim, std::vector& xVals, std::vector& temperatures); + ClassDefNV(Temperature, 1); }; @@ -270,8 +335,7 @@ struct Temperature { /// struct HV { - HV() - noexcept; + HV() noexcept; // Exmple strings // TPC_HV_A03_I_G1B_I @@ -483,5 +547,77 @@ struct Gas { ClassDefNV(Gas, 1); }; +struct RobustPressure { + using Stats = o2::math_utils::RollingStats; + Stats surfaceAtmosPressure; ///< rolling statistics of surface sensor + Stats cavernAtmosPressure; ///< rolling statistics of cavern sensor 1 + Stats cavernAtmosPressure2; ///< rolling statistics of cavern sensor 2 + Stats cavernAtmosPressure12; ///< rolling statistics of cavernAtmosPressure/cavernAtmosPressure2 + Stats cavernAtmosPressure1S; ///< rolling statistics of cavernAtmosPressure/surfaceAtmosPressure + Stats cavernAtmosPressure2S; ///< rolling statistics of cavernAtmosPressure2/surfaceAtmosPressure + std::vector isOk; ///< bit mask of valid sensors: cavernBit 0, cavern2Bit = 1, surfaceBit = 2 + std::vector robustPressure; ///< combined robust pressure value that should be used + std::vector time; ///< time stamps of all pressure values + TimeStampType timeInterval; ///< time interval used for rolling statistics + TimeStampType timeIntervalRef; ///< reference time interval used for normalization of pressure sensors + float maxDist{}; ///< maximum allowed time distance between sensors to be accepted for robust pressure calculation + float maxDiff{0.2f}; ///< maximum allowed pressure difference between sensors to be accepted for robust pressure calculation + + ClassDefNV(RobustPressure, 1); +}; + +struct Pressure { + + /// \brief fill pressure data + /// \param sensor name of the sensor from DCS data stream + /// \param time measurement time + /// \param value pressure value + void fill(std::string_view sensor, const TimeStampType time, const DataType value); + + /// sort pressure values and remove obvious outliers + /// \param pMin min accepted pressure + /// \param pMax max accepted pressure + void sortAndClean(float pMin = 800, float pMax = 1100); + + /// \clear all stored data except the buffer + void clear(); + + /// append other pressure values + void append(const Pressure& other); + + /// \return get minimum time of stored data + TimeStampType getMinTime() const; + + /// \return get maximum time of stored data + TimeStampType getMaxTime() const; + + /// \brief average pressure values for given time interval + /// \param timeInterval time interval for which the pressure values are averaged + /// \param timeIntervalRef time interval used to calculate the normalization values for the pressure + /// \param tStart min time of the data + /// \param tEnd max time of the data + /// \param nthreads numbe rof threads used for some calculations + void makeRobustPressure(TimeStampType timeInterval = 100 * 1000, TimeStampType timeIntervalRef = 24 * 60 * 1000, TimeStampType tStart = 1, TimeStampType tEnd = 0, const int nthreads = 1); + + /// set aliases for the cuts used in the calculation of the robust pressure + static void setAliases(TTree* tree); + + RawDPsF cavernAtmosPressure{}; ///< raw pressure in the cavern from sensor 1 + RawDPsF cavernAtmosPressure2{}; ///< raw pressure in the cavern from sensor 2 + RawDPsF surfaceAtmosPressure{}; ///< raw pressure at the surface + RobustPressure robustPressure{}; ///< combined robust pressure estimator from all three sensors + + std::pair, std::vector> mCavernAtmosPressure1Buff{}; ///, std::vector> mCavernAtmosPressure2Buff{}; ///, std::vector> mSurfaceAtmosPressureBuff{}; ///, std::vector> mPressure12Buff{}; ///, std::vector> mPressure1SBuff{}; ///, std::vector> mPressure2SBuff{}; ///, std::vector> mRobPressureBuff{}; /// #include "DataFormatsTPC/DCS.h" +#include "TLinearFitter.h" +#include "TTree.h" using namespace o2::tpc::dcs; @@ -117,3 +119,416 @@ TimeStampType Gas::getMaxTime() const return *std::max_element(times.begin(), times.end()); } + +TimeStampType Pressure::getMinTime() const +{ + constexpr auto max = std::numeric_limits::max(); + const std::vector times{ + cavernAtmosPressure.data.size() ? cavernAtmosPressure.data.front().time : max, + cavernAtmosPressure2.data.size() ? cavernAtmosPressure2.data.front().time : max, + surfaceAtmosPressure.data.size() ? surfaceAtmosPressure.data.front().time : max, + }; + + return *std::min_element(times.begin(), times.end()); +} + +TimeStampType Pressure::getMaxTime() const +{ + constexpr auto min = 0; + const std::vector times{ + cavernAtmosPressure.data.size() ? cavernAtmosPressure.data.back().time : min, + cavernAtmosPressure2.data.size() ? cavernAtmosPressure2.data.back().time : min, + surfaceAtmosPressure.data.size() ? surfaceAtmosPressure.data.back().time : min, + }; + + return *std::max_element(times.begin(), times.end()); +} + +bool Temperature::makeFit(TLinearFitter& fitter, const int nDim, std::vector& xVals, std::vector& temperatures) +{ + const int minPointsForFit = 5; + if (temperatures.empty() || (temperatures.size() < minPointsForFit)) { + LOGP(warning, "Number of points {} for fit smaller than minimum of {}!", temperatures.size(), minPointsForFit); + return false; + } + + fitter.ClearPoints(); + fitter.AssignData(temperatures.size(), nDim, xVals.data(), temperatures.data()); + int status = fitter.Eval(); + if (status == 1) { + LOGP(warning, "Fit failed!"); + return false; + } + return true; +} + +void Temperature::fitTemperature(Side side, dcs::TimeStampType fitInterval, const bool roundToInterval) +{ + // clear old data + auto& stats = (side == Side::A) ? statsA : statsC; + stats.clear(); + + // temperature fits in x-y + const int nDim = 2; + TLinearFitter fitter(nDim, "1 ++ x0 ++ x1", ""); + std::array startPos{}; + const size_t sensorOffset = (side == Side::C) ? dcs::Temperature::SensorsPerSide : 0; + + const dcs::TimeStampType refTime = getMinTime(raw, fitInterval, roundToInterval); + const dcs::TimeStampType refTimeMax = getMaxTime(raw); + + // calculate number of intervals and see if the last interval should be merged into the previous one + const int lastIntervalDuration = (refTimeMax - refTime) % fitInterval; + + // process the last interval only if it contains more than 50% of the interval duration + const bool procLastInt = (lastIntervalDuration / fitInterval > 0.5); + int numIntervals = (refTimeMax - refTime) / fitInterval + procLastInt; + if (numIntervals == 0) { + numIntervals = 1; + } + + // buffer for fit values + std::vector xVals; + std::vector temperatures; + xVals.reserve(2 * 1000); + temperatures.reserve(1000); + + for (int interval = 0; interval < numIntervals; ++interval) { + const dcs::TimeStampType timeStart = refTime + interval * fitInterval; + + // clear buffer + xVals.clear(); + temperatures.clear(); + + // TODO: check if we should use refTime + dcs::TimeStampType firstTime = std::numeric_limits::max(); + dcs::TimeStampType LastTime = 0; + + for (size_t iSensor = 0; iSensor < dcs::Temperature::SensorsPerSide; ++iSensor) { + const auto& sensor = raw[iSensor + sensorOffset]; + + LOGP(debug, "sensor {}, start {}, size {}", sensor.sensorNumber, startPos[iSensor], sensor.data.size()); + while (startPos[iSensor] < sensor.data.size()) { + const auto& dataPoint = sensor.data[startPos[iSensor]]; + if (((dataPoint.time - timeStart) >= fitInterval) && (interval != numIntervals - 1)) { + LOGP(debug, "sensor {}, {} - {} >= {}", sensor.sensorNumber, dataPoint.time, timeStart, fitInterval); + break; + } + firstTime = std::min(firstTime, dataPoint.time); + LastTime = std::max(LastTime, dataPoint.time); + const auto temperature = dataPoint.value; + // sanity check + ++startPos[iSensor]; + if (temperature < 15 || temperature > 25) { + continue; + } + const auto& pos = dcs::Temperature::SensorPosition[iSensor + sensorOffset]; + xVals.emplace_back(pos.x); + xVals.emplace_back(pos.y); + temperatures.emplace_back(temperature); + } + } + if (firstTime < std::numeric_limits::max() && !temperatures.empty()) { + const bool fitOk = makeFit(fitter, nDim, xVals, temperatures); + if (!fitOk) { + continue; + } + auto& stat = stats.data.emplace_back(); + stat.time = (firstTime + LastTime) / 2; + stat.value.mean = fitter.GetParameter(0); + stat.value.gradX = fitter.GetParameter(1); + stat.value.gradY = fitter.GetParameter(2); + + // check if data contains outliers + const float maxDeltaT = 1; + const float meanTemp = fitter.GetParameter(0); + const bool isDataGood = std::all_of(temperatures.begin(), temperatures.end(), [meanTemp, maxDeltaT](double t) { return std::abs(t - meanTemp) < maxDeltaT; }); + + // do second iteration only in case of outliers + if (!isDataGood) { + std::vector xVals2; + std::vector temperatures2; + xVals2.reserve(xVals.size()); + temperatures2.reserve(temperatures.size()); + for (int i = 0; i < temperatures.size(); ++i) { + if (std::abs(temperatures[i] - meanTemp) < maxDeltaT) { + const int idx = 2 * i; + xVals2.emplace_back(xVals[idx]); + xVals2.emplace_back(xVals[idx + 1]); + temperatures2.emplace_back(temperatures[i]); + } + } + const bool fitOk2 = makeFit(fitter, nDim, xVals2, temperatures2); + if (fitOk2) { + stat.value.mean = fitter.GetParameter(0); + stat.value.gradX = fitter.GetParameter(1); + stat.value.gradY = fitter.GetParameter(2); + } + } + } + } +} + +void Pressure::fill(std::string_view sensor, const TimeStampType time, const DataType value) +{ + if (sensor == "CavernAtmosPressure") { + cavernAtmosPressure.fill(time, value); + } else if (sensor == "CavernAtmosPressure2") { + cavernAtmosPressure2.fill(time, value); + } else if (sensor == "SurfaceAtmosPressure") { + surfaceAtmosPressure.fill(time, value); + } else { + LOGP(warning, "Unknown pressure sensor {}", sensor); + } +} + +void Pressure::sortAndClean(float pMin, float pMax) +{ + cavernAtmosPressure.sortAndClean(); + cavernAtmosPressure2.sortAndClean(); + surfaceAtmosPressure.sortAndClean(); + + auto removeOutliers = [](auto& dataVec, auto minVal, auto maxVal) { + dataVec.erase( + std::remove_if(dataVec.begin(), dataVec.end(), + [minVal, maxVal](const auto& dp) { + return (dp.value < minVal || dp.value > maxVal); + }), + dataVec.end()); + }; + + removeOutliers(cavernAtmosPressure.data, pMin, pMax); + removeOutliers(cavernAtmosPressure2.data, pMin, pMax); + removeOutliers(surfaceAtmosPressure.data, pMin, pMax); +} + +void Pressure::clear() +{ + cavernAtmosPressure.clear(); + cavernAtmosPressure2.clear(); + surfaceAtmosPressure.clear(); + robustPressure = RobustPressure(); +} + +void Pressure::append(const Pressure& other) +{ + cavernAtmosPressure.append(other.cavernAtmosPressure); + cavernAtmosPressure2.append(other.cavernAtmosPressure2); + surfaceAtmosPressure.append(other.surfaceAtmosPressure); +} + +void fillBuffer(std::pair, std::vector>& buffer, const std::pair, std::vector>& values, TimeStampType tStart, const int minPoints) +{ + const auto itStartBuff = std::lower_bound(buffer.second.begin(), buffer.second.end(), tStart); + size_t idxStartBuffer = std::distance(buffer.second.begin(), itStartBuff); + if (buffer.first.size() - idxStartBuffer < minPoints) { + if (buffer.first.size() < minPoints) { + idxStartBuffer = 0; + } else { + idxStartBuffer = buffer.first.size() - minPoints; + } + } + + std::pair, std::vector> buffTmp{ + std::vector(buffer.first.begin() + idxStartBuffer, buffer.first.end()), + std::vector(buffer.second.begin() + idxStartBuffer, buffer.second.end())}; + + buffTmp.first.insert(buffTmp.first.end(), values.first.begin(), values.first.end()); + buffTmp.second.insert(buffTmp.second.end(), values.second.begin(), values.second.end()); + + buffer = std::move(buffTmp); +} + +void Pressure::makeRobustPressure(TimeStampType timeInterval, TimeStampType timeIntervalRef, TimeStampType tStart, TimeStampType tEnd, const int nthreads) +{ + const auto surfaceAtmosPressurePair = surfaceAtmosPressure.getPairOfVector(); + const auto cavernAtmosPressurePair = cavernAtmosPressure.getPairOfVector(); + const auto cavernAtmosPressure2Pair = cavernAtmosPressure2.getPairOfVector(); + + // round to second + tStart = tStart / 1000 * 1000; + const TimeStampType tStartRef = (tStart - timeIntervalRef); + const int minPointsRef = 50; + fillBuffer(mCavernAtmosPressure1Buff, cavernAtmosPressurePair, tStartRef, minPointsRef); + fillBuffer(mCavernAtmosPressure2Buff, cavernAtmosPressure2Pair, tStartRef, minPointsRef); + fillBuffer(mSurfaceAtmosPressureBuff, surfaceAtmosPressurePair, tStartRef, minPointsRef); + + int nIntervals = std::round((tEnd - tStart) / timeInterval); + if (nIntervals == 0) { + nIntervals = 1; // at least one interval + } + std::vector times; + times.reserve(nIntervals); + for (int i = 0; i < nIntervals; ++i) { + times.emplace_back(tStart + (i + 0.5) * timeInterval); + } + + /// minimum number of points in the interval - otherwise use the n closest points + const int minPoints = 4; + const auto cavernAtmosPressureStats = o2::math_utils::getRollingStatistics(mCavernAtmosPressure1Buff.second, mCavernAtmosPressure1Buff.first, times, timeInterval, nthreads, minPoints, minPoints); + const auto cavernAtmosPressure2Stats = o2::math_utils::getRollingStatistics(mCavernAtmosPressure2Buff.second, mCavernAtmosPressure2Buff.first, times, timeInterval, nthreads, minPoints, minPoints); + const auto surfaceAtmosPressureStats = o2::math_utils::getRollingStatistics(mSurfaceAtmosPressureBuff.second, mSurfaceAtmosPressureBuff.first, times, timeInterval, nthreads, minPoints, minPoints); + + // subtract the moving median values from the different sensors if they are ok + std::pair, std::vector> cavernAtmosPressure12; + std::pair, std::vector> cavernAtmosPressure1S; + std::pair, std::vector> cavernAtmosPressure2S; + cavernAtmosPressure12.first.reserve(nIntervals); + cavernAtmosPressure1S.first.reserve(nIntervals); + cavernAtmosPressure2S.first.reserve(nIntervals); + cavernAtmosPressure12.second.reserve(nIntervals); + cavernAtmosPressure1S.second.reserve(nIntervals); + cavernAtmosPressure2S.second.reserve(nIntervals); + + for (int i = 0; i < nIntervals; i++) { + // coarse check if data is close by + const int maxDist = 600 * 1000; + const bool cavernOk = (cavernAtmosPressureStats.median[i] > 0) && (cavernAtmosPressureStats.closestDistanceL[i] < maxDist) && (cavernAtmosPressureStats.closestDistanceR[i] < maxDist); + const bool cavern2Ok = (cavernAtmosPressure2Stats.median[i] > 0) && (cavernAtmosPressure2Stats.closestDistanceL[i] < maxDist) && (cavernAtmosPressure2Stats.closestDistanceR[i] < maxDist); + const bool surfaceOk = (surfaceAtmosPressureStats.median[i] > 0) && (surfaceAtmosPressureStats.closestDistanceL[i] < maxDist) && (surfaceAtmosPressureStats.closestDistanceR[i] < maxDist); + + if (cavernOk && cavern2Ok) { + cavernAtmosPressure12.first.emplace_back(cavernAtmosPressureStats.median[i] - cavernAtmosPressure2Stats.median[i]); + cavernAtmosPressure12.second.emplace_back(times[i]); + } + if (cavernOk && surfaceOk) { + cavernAtmosPressure1S.first.emplace_back(cavernAtmosPressureStats.median[i] - surfaceAtmosPressureStats.median[i]); + cavernAtmosPressure1S.second.emplace_back(times[i]); + } + if (cavern2Ok && surfaceOk) { + cavernAtmosPressure2S.first.emplace_back(cavernAtmosPressure2Stats.median[i] - surfaceAtmosPressureStats.median[i]); + cavernAtmosPressure2S.second.emplace_back(times[i]); + } + } + + fillBuffer(mPressure12Buff, cavernAtmosPressure12, tStartRef, minPointsRef); + fillBuffer(mPressure1SBuff, cavernAtmosPressure1S, tStartRef, minPointsRef); + fillBuffer(mPressure2SBuff, cavernAtmosPressure2S, tStartRef, minPointsRef); + + // get long term median of diffs - this is used for normalization of the pressure values - + const auto cavernAtmosPressure12Stats = o2::math_utils::getRollingStatistics(mPressure12Buff.second, mPressure12Buff.first, times, timeIntervalRef, nthreads, 3, minPointsRef); + const auto cavernAtmosPressure1SStats = o2::math_utils::getRollingStatistics(mPressure1SBuff.second, mPressure1SBuff.first, times, timeIntervalRef, nthreads, 3, minPointsRef); + const auto cavernAtmosPressure2SStats = o2::math_utils::getRollingStatistics(mPressure2SBuff.second, mPressure2SBuff.first, times, timeIntervalRef, nthreads, 3, minPointsRef); + + // calculate diffs of median values + const float maxDist = 20 * timeInterval; + const float maxDiff = 0.2; + std::pair, std::vector> robustPressureTmp; + robustPressureTmp.first.reserve(nIntervals); + robustPressureTmp.second.reserve(nIntervals); + std::vector isOk(nIntervals); + + for (int i = 0; i < nIntervals; ++i) { + // difference beween pressure values corrected for the long term median + const float delta12 = cavernAtmosPressureStats.median[i] - cavernAtmosPressure2Stats.median[i] - cavernAtmosPressure12Stats.median[i]; + const float delta1S = cavernAtmosPressureStats.median[i] - surfaceAtmosPressureStats.median[i] - cavernAtmosPressure1SStats.median[i]; + const float delta2S = cavernAtmosPressure2Stats.median[i] - surfaceAtmosPressureStats.median[i] - cavernAtmosPressure2SStats.median[i]; + + const auto distCavernAtmosPressureL = cavernAtmosPressureStats.closestDistanceL[i]; + const auto distCavernAtmosPressure2L = cavernAtmosPressure2Stats.closestDistanceL[i]; + const auto distSurfaceAtmosPressureL = surfaceAtmosPressureStats.closestDistanceL[i]; + const auto distCavernAtmosPressureR = cavernAtmosPressureStats.closestDistanceR[i]; + const auto distCavernAtmosPressure2R = cavernAtmosPressure2Stats.closestDistanceR[i]; + const auto distSurfaceAtmosPressureR = surfaceAtmosPressureStats.closestDistanceR[i]; + + // check if data is ok + const bool cavernDistOk = (cavernAtmosPressureStats.median[i] > 0) && ((distCavernAtmosPressureL < maxDist) || (distCavernAtmosPressureR < maxDist)); + const bool cavern2DistOk = (cavernAtmosPressure2Stats.median[i] > 0) && ((distCavernAtmosPressure2L < maxDist) || (distCavernAtmosPressure2R < maxDist)); + const bool surfaceDistOk = (surfaceAtmosPressureStats.median[i] > 0) && ((distSurfaceAtmosPressureL < maxDist) || (distSurfaceAtmosPressureR < maxDist)); + const bool onlyOneSensor = (cavernDistOk + cavern2DistOk + surfaceDistOk) == 1; // check if only 1 sensor exists, if so use that sensor + + uint8_t maskIsOkTmp = 0; + const int cavernBit = 0; // val 1 + const int cavern2Bit = 1; // val 2 + const int surfaceBit = 2; // val 4 + + // check if ratio sensor 1 and 2 are good + // maskIsOkTmp = 3 + if (((std::abs(delta12) < maxDiff) && (cavernDistOk && cavern2DistOk)) || onlyOneSensor) { + if (cavernDistOk) { + maskIsOkTmp |= (1 << cavernBit); + } + if (cavern2DistOk) { + maskIsOkTmp |= (1 << cavern2Bit); + } + } + + // check if ratio sensor 1 and surface are good + // maskIsOkTmp = 5 + if ((std::abs(delta1S) < maxDiff) && ((cavernDistOk && surfaceDistOk)) || onlyOneSensor) { + if (cavernDistOk) { + maskIsOkTmp |= (1 << cavernBit); + } + if (surfaceDistOk) { + maskIsOkTmp |= (1 << surfaceBit); + } + } + + // check if ratio sensor 2 and surface are good + // maskIsOkTmp = 6 + if ((std::abs(delta2S) < maxDiff) && ((cavern2DistOk && surfaceDistOk)) || onlyOneSensor) { + if (cavern2DistOk) { + maskIsOkTmp |= (1 << cavern2Bit); + } + if (surfaceDistOk) { + maskIsOkTmp |= (1 << surfaceBit); + } + } + + // calculate robust pressure + float pressure = 0; + int pressureCount = 0; + if ((maskIsOkTmp >> cavernBit) & 1) { + pressure += cavernAtmosPressureStats.median[i]; + pressureCount++; + } + + if ((maskIsOkTmp >> cavern2Bit) & 1) { + pressure += cavernAtmosPressure2Stats.median[i] + cavernAtmosPressure12Stats.median[i]; + pressureCount++; + } + + if ((maskIsOkTmp >> surfaceBit) & 1) { + pressure += surfaceAtmosPressureStats.median[i] + cavernAtmosPressure1SStats.median[i]; + pressureCount++; + } + + isOk[i] = maskIsOkTmp; + if (pressureCount > 0) { + pressure /= pressureCount; + robustPressureTmp.first.emplace_back(pressure); + robustPressureTmp.second.emplace_back(times[i]); + } + } + + fillBuffer(mRobPressureBuff, robustPressureTmp, tStartRef, minPointsRef); + + RobustPressure& pOut = robustPressure; + pOut.surfaceAtmosPressure = std::move(surfaceAtmosPressureStats); + pOut.cavernAtmosPressure2 = std::move(cavernAtmosPressure2Stats); + pOut.cavernAtmosPressure = std::move(cavernAtmosPressureStats); + pOut.cavernAtmosPressure12 = std::move(cavernAtmosPressure12Stats); + pOut.cavernAtmosPressure1S = std::move(cavernAtmosPressure1SStats); + pOut.cavernAtmosPressure2S = std::move(cavernAtmosPressure2SStats); + pOut.isOk = std::move(isOk); + pOut.robustPressure = o2::math_utils::getRollingStatistics(mRobPressureBuff.second, mRobPressureBuff.first, times, timeInterval, nthreads, 1, 5).median; + pOut.time = std::move(times); + pOut.timeInterval = timeInterval; + pOut.timeIntervalRef = timeIntervalRef; + pOut.maxDist = maxDist; + pOut.maxDiff = maxDiff; +} + +void Pressure::setAliases(TTree* tree) +{ + tree->SetAlias("cavernDistOk", "robustPressure.cavernAtmosPressure.median>0 && (robustPressure.cavernAtmosPressure.closestDistanceRSetAlias("cavern2DistOk", "robustPressure.cavernAtmosPressure2.median>0 && (robustPressure.cavernAtmosPressure2.closestDistanceRSetAlias("surfaceDistOk", "robustPressure.surfaceAtmosPressure.median>0 && (robustPressure.surfaceAtmosPressure.closestDistanceRSetAlias("onlyOneSensor", "(cavernDistOk + cavern2DistOk + surfaceDistOk) == 1"); + tree->SetAlias("delta12", "robustPressure.cavernAtmosPressure.median - robustPressure.cavernAtmosPressure2.median - robustPressure.cavernAtmosPressure12.median"); + tree->SetAlias("delta1S", "robustPressure.cavernAtmosPressure.median - robustPressure.surfaceAtmosPressure.median - robustPressure.cavernAtmosPressure1S.median"); + tree->SetAlias("delta2S", "robustPressure.surfaceAtmosPressure.median - robustPressure.cavernAtmosPressure2.median - robustPressure.cavernAtmosPressure2S.median"); + tree->SetAlias("delta12_Ok", "abs(delta12)SetAlias("delta1S_Ok", "abs(delta1S)SetAlias("delta2S_Ok", "abs(delta2S) + ; #pragma link C++ class o2::tpc::dcs::DataPointVector < o2::tpc::dcs::HV::StackState> + ; #pragma link C++ class o2::tpc::dcs::Gas + ; +#pragma link C++ class o2::tpc::dcs::Pressure + ; +#pragma link C++ class o2::tpc::dcs::RobustPressure + ; #pragma link C++ class o2::tpc::PIDResponse + ; #pragma link C++ class o2::tpc::TriggerWordDLBZS + ; #pragma link C++ class o2::tpc::TriggerInfoDLBZS + ; diff --git a/Detectors/DCS/testWorkflow/src/DCSDataReplaySpec.cxx b/Detectors/DCS/testWorkflow/src/DCSDataReplaySpec.cxx index 783f6ae76e707..8dc003fc176f3 100644 --- a/Detectors/DCS/testWorkflow/src/DCSDataReplaySpec.cxx +++ b/Detectors/DCS/testWorkflow/src/DCSDataReplaySpec.cxx @@ -46,6 +46,10 @@ class DCSDataReplayer : public o2::framework::Task char mAlias[50]; uint64_t mMaxTF; uint64_t mTFs = 0; + int deltaTimeSendData = -1; + uint32_t startTime = -1; + uint32_t endTime = 0; + std::vector> dataIndicesPerTF; TTree mInputData; std::vector mDataPointHints; o2::header::DataDescription mDataDescription; @@ -59,6 +63,7 @@ void DCSDataReplayer::init(o2::framework::InitContext& ic) { mMaxTF = ic.options().get("max-timeframes"); mInputFileName = ic.options().get("input-file"); + deltaTimeSendData = ic.options().get("delta-time-send-data"); mInputData.ReadFile(mInputFileName.data(), "time/D:alias/C:value/D", ';'); mInputData.SetBranchAddress("time", &mTime); mInputData.SetBranchAddress("value", &mValue); @@ -73,16 +78,61 @@ void DCSDataReplayer::run(o2::framework::ProcessingContext& pc) LOG(info) << "Data generator reached TF " << tfid << ", stopping"; pc.services().get().endOfStream(); pc.services().get().readyToQuit(o2::framework::QuitRequest::Me); + return; } std::vector dpcoms; + for (Long64_t iEntry = 0; iEntry < mInputData.GetEntries(); ++iEntry) { - mInputData.GetEntry(iEntry); + int entryTree = iEntry; + + // load only releavant entries if requested + if (deltaTimeSendData > 0 && tfid > 2) { + + if (tfid - 1 >= dataIndicesPerTF.size()) { + LOGP(warning, "TF ID {} is larger than the number of TFs in dataIndicesPerTF: {}", tfid, dataIndicesPerTF.size()); + break; + } + + if (iEntry >= dataIndicesPerTF[tfid - 1].size()) { + break; + } else { + entryTree = dataIndicesPerTF[tfid - 1][iEntry]; + } + } + + mInputData.GetEntry(entryTree); const auto ultime = uint64_t(std::round(mTime * 1000)); const auto seconds = uint32_t(ultime / 1000); const auto msec = uint16_t(ultime % 1000); - - dpcoms.emplace_back(o2::dcs::createDataPointCompositeObject(mAlias, float(mValue), seconds, msec)); + if (deltaTimeSendData > 0) { + // send data in packages + if (tfid == 0) { + startTime = std::min(startTime, seconds); + endTime = std::max(endTime, seconds); + if (iEntry == mInputData.GetEntries() - 1) { + const int totalTFs = (endTime - startTime) / deltaTimeSendData + 1; + dataIndicesPerTF.resize(totalTFs); + LOGP(info, "Sending data from {} to {} with {} TFs", startTime, endTime, totalTFs); + } + } else { + if (tfid == 1) { + const int index = (seconds - startTime) / deltaTimeSendData; + dataIndicesPerTF[index].emplace_back(iEntry); + } + const uint64_t startTimeTF = startTime + (tfid - 1) * deltaTimeSendData; + const uint64_t endTimeTF = startTimeTF + deltaTimeSendData; + if (seconds >= startTimeTF && seconds < endTimeTF) { + dpcoms.emplace_back(o2::dcs::createDataPointCompositeObject(mAlias, float(mValue), seconds, msec)); + // check if all data has been processed + if (seconds == endTime) { + mMaxTF = tfid; + } + } + } + } else { + dpcoms.emplace_back(o2::dcs::createDataPointCompositeObject(mAlias, float(mValue), seconds, msec)); + } } // auto dpcoms = generate(mDataPointHints, fraction, tfid); @@ -113,6 +163,7 @@ o2::framework::DataProcessorSpec getDCSDataReplaySpec(std::vector CDBTypeMap{ {CDBType::CalGas, "TPC/Calib/Gas"}, {CDBType::CalTemperature, "TPC/Calib/Temperature"}, {CDBType::CalHV, "TPC/Calib/HV"}, + {CDBType::CalPressure, "TPC/Calib/Pressure"}, {CDBType::CalTopologyGain, "TPC/Calib/TopologyGainPiecewise"}, {CDBType::CalVDriftTgl, "TPC/Calib/VDriftTgl"}, // diff --git a/Detectors/TPC/calibration/CMakeLists.txt b/Detectors/TPC/calibration/CMakeLists.txt index 7722fc4e2884f..8bcb3254edb32 100644 --- a/Detectors/TPC/calibration/CMakeLists.txt +++ b/Detectors/TPC/calibration/CMakeLists.txt @@ -58,6 +58,7 @@ o2_add_library(TPCCalibration src/TPCMShapeCorrection.cxx src/DigitAdd.cxx src/CorrectdEdxDistortions.cxx + src/PressureTemperatureHelper.cxx PUBLIC_LINK_LIBRARIES O2::DataFormatsTPC O2::TPCBase O2::TPCReconstruction ROOT::Minuit Microsoft.GSL::GSL @@ -115,7 +116,8 @@ o2_target_root_dictionary(TPCCalibration include/TPCCalibration/CorrMapParam.h include/TPCCalibration/TPCMShapeCorrection.h include/TPCCalibration/DigitAdd.h - include/TPCCalibration/CorrectdEdxDistortions.h) + include/TPCCalibration/CorrectdEdxDistortions.h + include/TPCCalibration/PressureTemperatureHelper.h) o2_add_test_root_macro(macro/comparePedestalsAndNoise.C PUBLIC_LINK_LIBRARIES O2::TPCBase diff --git a/Detectors/TPC/calibration/include/TPCCalibration/PressureTemperatureHelper.h b/Detectors/TPC/calibration/include/TPCCalibration/PressureTemperatureHelper.h new file mode 100644 index 0000000000000..b636fdd2f296d --- /dev/null +++ b/Detectors/TPC/calibration/include/TPCCalibration/PressureTemperatureHelper.h @@ -0,0 +1,84 @@ +// 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 PressureTemperatureHelper.h +/// \brief Helper class to extract pressure and temperature +/// \author Matthias Kleiner + +#ifndef PRESSURETEMPERATUREHELPER_H_ +#define PRESSURETEMPERATUREHELPER_H_ + +#include "GPUCommonRtypes.h" +#include "Headers/DataHeader.h" +#include "CommonDataFormat/Pair.h" + +namespace o2::framework +{ +class ProcessingContext; +class ConcreteDataMatcher; +class InputSpec; +class OutputSpec; +} // namespace o2::framework + +namespace o2::tpc +{ + +class PressureTemperatureHelper +{ + public: + PressureTemperatureHelper() = default; + + /// check for new CCDB objects + bool accountCCDBInputs(const o2::framework::ConcreteDataMatcher& matcher, void* obj); + + /// trigger checking for CCDB objects + void extractCCDBInputs(o2::framework::ProcessingContext& pc) const; + + // add required inputs + static void requestCCDBInputs(std::vector& inputs); + + /// define outputs in case pressure and temperature will be send + static void setOutputs(std::vector& outputs); + + /// send temperature and pressure for given time stamp + void sendPTForTS(o2::framework::ProcessingContext& pc, const uint64_t timestamp) const; + + /// set fit interval range for temperature in ms + void setFitIntervalTemp(const int fitIntervalMS) { mFitIntervalMS = fitIntervalMS; } + + /// \brief interpolate input values for given timestamp + /// \param timestamps time stamps of the data + /// \param values data points + /// \param timestamp time where to interpolate the values + float interpolate(const std::vector& timestamps, const std::vector& values, uint64_t timestamp) const; + + /// get pressure for given time stamp in ms + float getPressure(const uint64_t timestamp) const { return interpolate(mPressure.second, mPressure.first, timestamp); } + + /// get temperature for given time stamp in ms + dataformats::Pair getTemperature(const uint64_t timestamp) const { return dataformats::Pair{interpolate(mTemperatureA.second, mTemperatureA.first, timestamp), interpolate(mTemperatureC.second, mTemperatureC.first, timestamp)}; } + + static constexpr o2::header::DataDescription getDataDescriptionPressure() { return o2::header::DataDescription{"pressure"}; } + static constexpr o2::header::DataDescription getDataDescriptionTemperature() { return o2::header::DataDescription{"temperature"}; } + + protected: + static void addInput(std::vector& inputs, o2::framework::InputSpec&& isp); + static void addOutput(std::vector& outputs, o2::framework::OutputSpec&& osp); + + std::pair, std::vector> mPressure; ///< pressure values for both measurements + std::pair, std::vector> mTemperatureA; ///< temperature values A-side + std::pair, std::vector> mTemperatureC; ///< temperature values C-side + int mFitIntervalMS{5 * 60 * 1000}; ///< fit interval for the temperature + + ClassDefNV(PressureTemperatureHelper, 1); +}; +} // namespace o2::tpc +#endif diff --git a/Detectors/TPC/calibration/src/PressureTemperatureHelper.cxx b/Detectors/TPC/calibration/src/PressureTemperatureHelper.cxx new file mode 100644 index 0000000000000..54145f0ecfaf1 --- /dev/null +++ b/Detectors/TPC/calibration/src/PressureTemperatureHelper.cxx @@ -0,0 +1,119 @@ +// 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 PressureTemperatureHelper.cxx +/// \brief Helper class to extract pressure and temperature +/// \author Matthias Kleiner + +#include "TPCCalibration/PressureTemperatureHelper.h" +#include "TPCBase/CDBInterface.h" +#include "Framework/ProcessingContext.h" +#include "DataFormatsTPC/DCS.h" +#include "Framework/InputRecord.h" +#include "Framework/CCDBParamSpec.h" +#include "Framework/DataAllocator.h" + +using namespace o2::tpc; +using namespace o2::framework; + +void PressureTemperatureHelper::extractCCDBInputs(ProcessingContext& pc) const +{ + pc.inputs().get("pressure"); + pc.inputs().get("temperature"); +} + +bool PressureTemperatureHelper::accountCCDBInputs(const ConcreteDataMatcher& matcher, void* obj) +{ + if (matcher == ConcreteDataMatcher(o2::header::gDataOriginTPC, "PRESSURECCDB", 0)) { + LOGP(info, "Updating pressure"); + const auto& pressure = ((dcs::Pressure*)obj); + mPressure.second = pressure->robustPressure.time; + mPressure.first = pressure->robustPressure.robustPressure; + return true; + } + + if (matcher == ConcreteDataMatcher(o2::header::gDataOriginTPC, "TEMPERATURECCDB", 0)) { + LOGP(info, "Updating temperature"); + auto temp = *(dcs::Temperature*)obj; + temp.fitTemperature(o2::tpc::Side::A, mFitIntervalMS, false); + temp.fitTemperature(o2::tpc::Side::C, mFitIntervalMS, false); + + mTemperatureA.first.clear(); + mTemperatureC.first.clear(); + mTemperatureA.second.clear(); + mTemperatureC.second.clear(); + + for (const auto& dp : temp.statsA.data) { + mTemperatureA.first.emplace_back(dp.value.mean); + mTemperatureA.second.emplace_back(dp.time); + } + + for (const auto& dp : temp.statsC.data) { + mTemperatureC.first.emplace_back(dp.value.mean); + mTemperatureC.second.emplace_back(dp.time); + } + return true; + } + return false; +} + +void PressureTemperatureHelper::requestCCDBInputs(std::vector& inputs) +{ + addInput(inputs, {"pressure", o2::header::gDataOriginTPC, "PRESSURECCDB", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CalPressure), {}, 1)}); + addInput(inputs, {"temperature", o2::header::gDataOriginTPC, "TEMPERATURECCDB", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CalTemperature), {}, 1)}); +} + +void PressureTemperatureHelper::addInput(std::vector& inputs, InputSpec&& isp) +{ + if (std::find(inputs.begin(), inputs.end(), isp) == inputs.end()) { + inputs.emplace_back(isp); + } +} + +void PressureTemperatureHelper::setOutputs(std::vector& outputs) +{ + addOutput(outputs, {o2::header::gDataOriginTPC, o2::tpc::PressureTemperatureHelper::getDataDescriptionPressure(), 0, Lifetime::Timeframe}); + addOutput(outputs, {o2::header::gDataOriginTPC, o2::tpc::PressureTemperatureHelper::getDataDescriptionTemperature(), 0, Lifetime::Timeframe}); +} + +void PressureTemperatureHelper::addOutput(std::vector& outputs, OutputSpec&& osp) +{ + if (std::find(outputs.begin(), outputs.end(), osp) == outputs.end()) { + outputs.emplace_back(osp); + } +} + +float PressureTemperatureHelper::interpolate(const std::vector& timestamps, const std::vector& values, uint64_t timestamp) const +{ + if (auto idxClosest = o2::math_utils::findClosestIndices(timestamps, timestamp)) { + auto [idxLeft, idxRight] = *idxClosest; + if (idxRight > idxLeft) { + const uint64_t x0 = timestamps[idxLeft]; + const uint64_t x1 = timestamps[idxRight]; + const float y0 = values[idxLeft]; + const float y1 = values[idxRight]; + const float y = (y0 * (x1 - timestamp) + y1 * (timestamp - x0)) / (x1 - x0); + return y; + } else { + return values[idxLeft]; + } + } + return 0; // this should never happen +} + +void PressureTemperatureHelper::sendPTForTS(o2::framework::ProcessingContext& pc, const uint64_t timestamp) const +{ + const float pressure = getPressure(timestamp); + const auto temp = getTemperature(timestamp); + LOGP(info, "Sending pressure {}, temperature A {} and temperature C {} for timestamp {}", pressure, temp.first, temp.second, timestamp); + pc.outputs().snapshot(Output{o2::header::gDataOriginTPC, o2::tpc::PressureTemperatureHelper::getDataDescriptionTemperature()}, temp); + pc.outputs().snapshot(Output{o2::header::gDataOriginTPC, o2::tpc::PressureTemperatureHelper::getDataDescriptionPressure()}, pressure); +} diff --git a/Detectors/TPC/calibration/src/TPCCalibrationLinkDef.h b/Detectors/TPC/calibration/src/TPCCalibrationLinkDef.h index d42627197cd7f..6e15e2dd0427a 100644 --- a/Detectors/TPC/calibration/src/TPCCalibrationLinkDef.h +++ b/Detectors/TPC/calibration/src/TPCCalibrationLinkDef.h @@ -122,4 +122,5 @@ #pragma link C++ struct o2::tpc::BoundaryPotentialIFC + ; #pragma link C++ class o2::tpc::DigitAdd + ; #pragma link C++ class std::vector < o2::tpc::DigitAdd> + ; +#pragma link C++ class o2::tpc::PressureTemperatureHelper + ; #endif diff --git a/Detectors/TPC/dcs/CMakeLists.txt b/Detectors/TPC/dcs/CMakeLists.txt index 9db4e26a7429d..31524dd5f2c2f 100644 --- a/Detectors/TPC/dcs/CMakeLists.txt +++ b/Detectors/TPC/dcs/CMakeLists.txt @@ -18,8 +18,7 @@ o2_add_library(TPCdcs PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsDCS O2::DataFormatsTPC - O2::TPCBase - ROOT::Minuit) + O2::TPCBase) o2_target_root_dictionary(TPCdcs HEADERS include/TPCdcs/DCSProcessor.h) diff --git a/Detectors/TPC/dcs/include/TPCdcs/DCSProcessor.h b/Detectors/TPC/dcs/include/TPCdcs/DCSProcessor.h index 23c6df26c0fd7..e6ead9b0cb302 100644 --- a/Detectors/TPC/dcs/include/TPCdcs/DCSProcessor.h +++ b/Detectors/TPC/dcs/include/TPCdcs/DCSProcessor.h @@ -46,22 +46,14 @@ class DCSProcessor void fillTemperature(const DPCOM& dp); void fillHV(const DPCOM& dp); void fillGas(const DPCOM& dp); + void fillPressure(const DPCOM& dp); void finalizeSlot(); void finalize(); void finalizeTemperature(); void finalizeHighVoltage(); void finalizeGas(); - - void fitTemperature(Side side); - - /// get minimum time over all sensors. Assumes data is sorted in time - template - dcs::TimeStampType getMinTime(const std::vector>& data); - - /// get maximum time over all sensors. Assumes data is sorted in time - template - dcs::TimeStampType getMaxTime(const std::vector>& data); + void finalizePressure(); /// name of the debug output tree void setDebugOutputName(std::string_view name) { mDebugOutputName = name; } @@ -75,9 +67,17 @@ class DCSProcessor /// set the fit interval void setFitInterval(dcs::TimeStampType interval) { mFitInterval = interval; } + /// set the interval for averaging the pressure values + void setPressureInterval(dcs::TimeStampType interval) { mPressureInterval = interval; } + + void setRefPressureInterval(dcs::TimeStampType interval) { mPressureIntervalRef = interval; } + /// get fit interval auto getFitInterval() const { return mFitInterval; } + /// get fit interval + auto getPressureInterval() const { return mPressureInterval; } + /// round to fit interval void setRoundToInterval(const bool round = true) { mRoundToInterval = round; } @@ -87,6 +87,7 @@ class DCSProcessor mTemperature.clear(); mHighVoltage.clear(); mGas.clear(); + mPressure.clear(); mTimeTemperature = {}; mTimeHighVoltage = {}; @@ -99,21 +100,27 @@ class DCSProcessor const auto& getTimeTemperature() const { return mTimeTemperature; } const auto& getTimeHighVoltage() const { return mTimeHighVoltage; } const auto& getTimeGas() const { return mTimeGas; } + const auto& getTimePressure() const { return mTimePressure; } auto& getTemperature() { return mTemperature; } auto& getHighVoltage() { return mHighVoltage; } auto& getGas() { return mGas; } + auto& getPressure() { return mPressure; } private: dcs::Temperature mTemperature; ///< temperature value store dcs::HV mHighVoltage; ///< HV value store dcs::Gas mGas; ///< Gas value store + dcs::Pressure mPressure; ///< Pressure value TimeRange mTimeTemperature; ///< Time range for temperature values TimeRange mTimeHighVoltage; ///< Time range for high voltage values TimeRange mTimeGas; ///< Time range for gas values + TimeRange mTimePressure; ///< Time range for pressure values dcs::TimeStampType mFitInterval{5 * 60 * 1000}; ///< fit interval (ms) e.g. for temparature data + dcs::TimeStampType mPressureInterval{200 * 1000}; ///< interval (ms) for averaging pressure values + dcs::TimeStampType mPressureIntervalRef{60 * 60 * 1000}; ///< interval (ms) for averaging pressure values for longer reference time interval bool mWriteDebug{false}; ///< switch to dump debug tree bool mRoundToInterval{false}; ///< round to full fit interval e.g. full minute bool mHasData{false}; ///< if there are data to process @@ -123,44 +130,5 @@ class DCSProcessor ClassDefNV(DCSProcessor, 0); }; -template -dcs::TimeStampType DCSProcessor::getMinTime(const std::vector>& data) -{ - constexpr auto max = std::numeric_limits::max(); - dcs::TimeStampType firstTime = std::numeric_limits::max(); - for (const auto& sensor : data) { - const auto time = sensor.data.size() ? sensor.data.front().time : max; - firstTime = std::min(firstTime, time); - } - - // mFitInterval is is seconds. Round to full amount. - // if e.g. mFitInterval = 5min, then round 10:07:20.510 to 10:05:00.000 - if (mRoundToInterval) { - firstTime -= (firstTime % mFitInterval); - } - - return firstTime; -} - -template -dcs::TimeStampType DCSProcessor::getMaxTime(const std::vector>& data) -{ - constexpr auto min = 0; - dcs::TimeStampType lastTime = 0; - for (const auto& sensor : data) { - const auto time = sensor.data.size() ? sensor.data.back().time : 0; - lastTime = std::max(lastTime, time); - } - - // mFitInterval is is seconds. Round to full amount. - // if e.g. mFitInterval = 5min, then round 10:07:20.510 to 10:05:00.000 - // TODO: fix this - // if (mRoundToInterval) { - // lastTime -= (lastTime % mFitInterval); - //} - - return lastTime; -} - } // namespace o2::tpc #endif diff --git a/Detectors/TPC/dcs/src/DCSDPHints.cxx b/Detectors/TPC/dcs/src/DCSDPHints.cxx index 00ca097a455eb..02ada7d588f9b 100644 --- a/Detectors/TPC/dcs/src/DCSDPHints.cxx +++ b/Detectors/TPC/dcs/src/DCSDPHints.cxx @@ -98,5 +98,10 @@ std::vector o2::tpc::dcs::getTPCDCSDPHints(const int ma dphints.emplace_back(o2::dcs::test::DataPointHint{fmt::format("TPC_HV_C[00..{:02}]_I_STATUS", maxSectors), 0, 29}); dphints.emplace_back(o2::dcs::test::DataPointHint{fmt::format("TPC_HV_C[00..{:02}]_O[1..3]_STATUS", maxSectors), 0, 29}); + // ===| pressure values |======================================================= + dphints.emplace_back(o2::dcs::test::DataPointHint{"CavernAtmosPressure", 850., 1050.}); + dphints.emplace_back(o2::dcs::test::DataPointHint{"CavernAtmosPressure2", 850., 1050.}); + dphints.emplace_back(o2::dcs::test::DataPointHint{"SurfaceAtmosPressure", 850., 1050.}); + return dphints; } diff --git a/Detectors/TPC/dcs/src/DCSProcessor.cxx b/Detectors/TPC/dcs/src/DCSProcessor.cxx index a26bce43e5c2e..4c9d196432687 100644 --- a/Detectors/TPC/dcs/src/DCSProcessor.cxx +++ b/Detectors/TPC/dcs/src/DCSProcessor.cxx @@ -15,10 +15,6 @@ #include -// ROOT includes -#include "TLinearFitter.h" -#include "TVectorD.h" - // O2 includes #include "DetectorsDCS/DataPointIdentifier.h" #include "DetectorsDCS/DataPointValue.h" @@ -43,6 +39,8 @@ void DCSProcessor::process(const gsl::span dps) constexpr auto HV_ID{"TPC_HV"sv}; constexpr auto GAS_ID1{"TPC_GC"sv}; constexpr auto GAS_ID2{"TPC_An"sv}; + constexpr auto PRESS_ID1{"Cavern"sv}; + constexpr auto PRESS_ID2{"Surfac"sv}; for (const auto& dp : dps) { const std::string_view alias(dp.id.get_alias()); @@ -56,8 +54,11 @@ void DCSProcessor::process(const gsl::span dps) } else if (id == GAS_ID1 || id == GAS_ID2) { LOGP(debug, "Gas DP: {}", alias); fillGas(dp); + } else if (id == PRESS_ID1 || id == PRESS_ID2) { + LOGP(debug, "Pressure DP: {}", alias); + fillPressure(dp); } else { - LOGP(warning, "Unknown data point: {}", alias); + LOGP(warning, "Unknown data point: {} with id {}", alias, id); } } } @@ -125,84 +126,38 @@ void DCSProcessor::fillGas(const DPCOM& dp) mGas.fill(alias, time, value); } +void DCSProcessor::fillPressure(const DPCOM& dp) +{ + const std::string_view alias(dp.id.get_alias()); + const auto value = getValueF(dp); + const auto time = dp.data.get_epoch_time(); + mPressure.fill(alias, time, value); +} + void DCSProcessor::finalizeSlot() { finalizeTemperature(); finalizeHighVoltage(); finalizeGas(); + finalizePressure(); mHasData = false; } -void DCSProcessor::fitTemperature(Side side) -{ - //// temperature fits in x-y - TLinearFitter fitter(3, "x0 ++ x1 ++ x2"); - bool nextInterval = true; - std::array startPos{}; - const size_t sensorOffset = (side == Side::C) ? dcs::Temperature::SensorsPerSide : 0; - dcs::TimeStampType refTime = getMinTime(mTemperature.raw); - - while (nextInterval) { - // TODO: check if we should use refTime - dcs::TimeStampType firstTime = std::numeric_limits::max(); - - nextInterval = false; - for (size_t iSensor = 0; iSensor < dcs::Temperature::SensorsPerSide; ++iSensor) { - const auto& sensor = mTemperature.raw[iSensor + sensorOffset]; - - LOGP(debug, "sensor {}, start {}, size {}", sensor.sensorNumber, startPos[iSensor], sensor.data.size()); - while (startPos[iSensor] < sensor.data.size()) { - const auto& dataPoint = sensor.data[startPos[iSensor]]; - if ((dataPoint.time - refTime) >= mFitInterval) { - LOGP(debug, "sensor {}, {} - {} >= {}", sensor.sensorNumber, dataPoint.time, refTime, mFitInterval); - break; - } - nextInterval = true; - firstTime = std::min(firstTime, dataPoint.time); - const auto temperature = dataPoint.value; - // sanity check - if (temperature < 15 || temperature > 25) { - ++startPos[iSensor]; - continue; - } - const auto& pos = dcs::Temperature::SensorPosition[iSensor + sensorOffset]; - double x[] = {1., double(pos.x), double(pos.y)}; - fitter.AddPoint(x, temperature, 1); - ++startPos[iSensor]; - } - } - if (firstTime < std::numeric_limits::max()) { - fitter.Eval(); - LOGP(info, "Side {}, fit interval {} - {} with {} points", int(side), refTime, refTime + mFitInterval - 1, fitter.GetNpoints()); - - auto& stats = (side == Side::A) ? mTemperature.statsA : mTemperature.statsC; - auto& stat = stats.data.emplace_back(); - stat.time = firstTime; - stat.value.mean = fitter.GetParameter(0); - stat.value.gradX = fitter.GetParameter(1); - stat.value.gradY = fitter.GetParameter(2); - - fitter.ClearPoints(); - refTime += mFitInterval; - } - } -} - void DCSProcessor::finalizeTemperature() { mTemperature.sortAndClean(); - fitTemperature(Side::A); - fitTemperature(Side::C); - mTimeTemperature = {getMinTime(mTemperature.raw), getMaxTime(mTemperature.raw)}; + mTemperature.fitTemperature(Side::A, mFitInterval, mRoundToInterval); + mTemperature.fitTemperature(Side::C, mFitInterval, mRoundToInterval); + mTimeTemperature = {getMinTime(mTemperature.raw, mRoundToInterval, mFitInterval), getMaxTime(mTemperature.raw)}; } void DCSProcessor::finalizeHighVoltage() { mHighVoltage.sortAndClean(); - auto minTime = getMinTime(mHighVoltage.currents); - minTime = std::min(minTime, getMinTime(mHighVoltage.voltages)); - minTime = std::min(minTime, getMinTime(mHighVoltage.states)); + auto minTime = getMinTime(mHighVoltage.currents, mRoundToInterval, mFitInterval); + minTime = std::min(minTime, getMinTime(mHighVoltage.voltages, mRoundToInterval, mFitInterval)); + minTime = std::min(minTime, getMinTime(mHighVoltage.states, mRoundToInterval, mFitInterval)); auto maxTime = getMaxTime(mHighVoltage.currents); maxTime = std::max(maxTime, getMaxTime(mHighVoltage.voltages)); @@ -217,6 +172,16 @@ void DCSProcessor::finalizeGas() mTimeGas = {mGas.getMinTime(), mGas.getMaxTime()}; } +void DCSProcessor::finalizePressure() +{ + mPressure.sortAndClean(); + mTimePressure = {mPressure.getMinTime(), mPressure.getMaxTime()}; + // if there is data perform the processing + if (mTimePressure.last > 0) { + mPressure.makeRobustPressure(mPressureInterval, mPressureIntervalRef, mTimePressure.first, mTimePressure.last); + } +} + void DCSProcessor::writeDebug() { if (!mDebugStream) { @@ -227,6 +192,7 @@ void DCSProcessor::writeDebug() << "Temperature=" << mTemperature << "HV=" << mHighVoltage << "Gas=" << mGas + << "Pressure=" << mPressure << "\n"; } diff --git a/Detectors/TPC/dcs/src/DCSSpec.cxx b/Detectors/TPC/dcs/src/DCSSpec.cxx index f99ff8f8aaaab..1b64ff7a75ba4 100644 --- a/Detectors/TPC/dcs/src/DCSSpec.cxx +++ b/Detectors/TPC/dcs/src/DCSSpec.cxx @@ -48,7 +48,7 @@ const std::unordered_map CDBDescMap{ {CDBType::CalTemperature, o2::header::DataDescription{"TPC_Temperature"}}, {CDBType::CalHV, o2::header::DataDescription{"TPC_HighVoltage"}}, {CDBType::CalGas, o2::header::DataDescription{"TPC_Gas"}}, -}; + {CDBType::CalPressure, o2::header::DataDescription{"TPC_Pressure"}}}; class DCSDevice : public o2::framework::Task { @@ -105,6 +105,7 @@ class DCSDevice : public o2::framework::Task bool mDebugWritten{false}; bool mWriteDebug{false}; bool mReportTiming{false}; + int mUpdateIntervalnTFs{-1}; }; void DCSDevice::init(o2::framework::InitContext& ic) @@ -112,6 +113,8 @@ void DCSDevice::init(o2::framework::InitContext& ic) mWriteDebug = ic.options().get("write-debug"); mCCDBupdateInterval = ic.options().get("update-interval"); mFitInterval = ic.options().get("fit-interval"); + const int pressureInterval = ic.options().get("pressure-interval"); + const int pressureIntervalRef = ic.options().get("pressure-ref-interval"); if (mCCDBupdateInterval < 0) { mCCDBupdateInterval = 0; } @@ -119,15 +122,18 @@ void DCSDevice::init(o2::framework::InitContext& ic) LOGP(info, "fit interval {} >= ccdb update interval {}, making them identical", mFitInterval, mCCDBupdateInterval); mFitInterval = mCCDBupdateInterval; } + mUpdateIntervalnTFs = ic.options().get("update-interval-nTFs"); mDCS.setFitInterval(mFitInterval * 1000); // in ms in mDCS + mDCS.setPressureInterval(pressureInterval * 1000); + mDCS.setRefPressureInterval(pressureIntervalRef * 1000); mDCS.setRoundToInterval(ic.options().get("round-to-interval")); // set default meta data mCDBStorage.setResponsible("Jens Wiechula (jens.wiechula@cern.ch)"); mCDBStorage.setIntervention(CDBIntervention::Automatic); mCDBStorage.setReason("DCS workflow upload"); - mReportTiming = ic.options().get("report-timing") || mWriteDebug; + mReportTiming = ic.options().get("report-timing"); } void DCSDevice::run(o2::framework::ProcessingContext& pc) @@ -137,10 +143,16 @@ void DCSDevice::run(o2::framework::ProcessingContext& pc) if (mUpdateIntervalStart == 0) { mUpdateIntervalStart = mLastCreationTime; } - if (mLastCreationTime - mUpdateIntervalStart >= uint64_t(mCCDBupdateInterval * 1000)) { + if (mUpdateIntervalnTFs > 0 && (pc.services().get().tfCounter % mUpdateIntervalnTFs == 0)) { + // finalize DCS for every n-TFs (useful for testing purpose when reading in data from local file) finalizeDCS(pc.outputs()); - mUpdateIntervalStart = mLastCreationTime; + } else { + if (mLastCreationTime - mUpdateIntervalStart >= uint64_t(mCCDBupdateInterval * 1000)) { + finalizeDCS(pc.outputs()); + mUpdateIntervalStart = mLastCreationTime; + } } + auto dps = pc.inputs().get>("input"); mDCS.process(dps); sw.Stop(); @@ -170,6 +182,7 @@ void DCSDevice::updateCCDB(DataAllocator& output) sendObject(output, mDCS.getTemperature(), CDBType::CalTemperature); sendObject(output, mDCS.getHighVoltage(), CDBType::CalHV); sendObject(output, mDCS.getGas(), CDBType::CalGas); + sendObject(output, mDCS.getPressure(), CDBType::CalPressure); } /// ===| create DCS processor |================================================= @@ -188,6 +201,9 @@ DataProcessorSpec getDCSSpec() outputs.emplace_back(ConcreteDataTypeMatcher{CDBPayload, CDBDescMap.at(CDBType::CalGas)}, Lifetime::Sporadic); outputs.emplace_back(ConcreteDataTypeMatcher{CDBWrapper, CDBDescMap.at(CDBType::CalGas)}, Lifetime::Sporadic); + outputs.emplace_back(ConcreteDataTypeMatcher{CDBPayload, CDBDescMap.at(CDBType::CalPressure)}, Lifetime::Sporadic); + outputs.emplace_back(ConcreteDataTypeMatcher{CDBWrapper, CDBDescMap.at(CDBType::CalPressure)}, Lifetime::Sporadic); + return DataProcessorSpec{ "tpc-dcs", Inputs{{"input", "DCS", "TPCDATAPOINTS"}}, @@ -199,7 +215,9 @@ DataProcessorSpec getDCSSpec() {"update-interval", VariantType::Int, 60 * 15, {"update interval in seconds for which ccdb entries are written"}}, {"fit-interval", VariantType::Int, 60 * 5, {"interval in seconds for which to e.g. perform fits of the temperature sensors"}}, {"round-to-interval", VariantType::Bool, false, {"round fit interval to fixed times e.g. to every 5min in the hour"}}, - } // end Options + {"pressure-interval", VariantType::Int, 100, {"interval in seconds for which to average the pressure values"}}, + {"pressure-ref-interval", VariantType::Int, 24 * 60 * 60, {"interval in seconds for which to calculate the reference pressure values"}}, + {"update-interval-nTFs", VariantType::Int, -1, {"only used when larger than 0: update interval in nTFs for which ccdb entries are written "}}} // end Options }; // end DataProcessorSpec } diff --git a/Detectors/TPC/workflow/CMakeLists.txt b/Detectors/TPC/workflow/CMakeLists.txt index 3b05e5067108c..fe7c9175968b5 100644 --- a/Detectors/TPC/workflow/CMakeLists.txt +++ b/Detectors/TPC/workflow/CMakeLists.txt @@ -44,6 +44,7 @@ o2_add_library(TPCWorkflow src/TPCTimeSeriesSpec.cxx src/TPCTimeSeriesWriterSpec.cxx src/TPCScalerSpec.cxx + src/TPCPressureTemperatureSpec.cxx TARGETVARNAME targetName PUBLIC_LINK_LIBRARIES O2::Framework O2::DataFormatsTPC O2::DPLUtils O2::TPCReconstruction @@ -293,4 +294,9 @@ o2_add_executable(refitter COMPONENT_NAME tpc PUBLIC_LINK_LIBRARIES O2::TPCWorkflowStudies) +o2_add_executable(pressure-temperature + COMPONENT_NAME tpc + SOURCES src/tpc-pressure-temperature.cxx + PUBLIC_LINK_LIBRARIES O2::TPCWorkflow) + add_subdirectory(readers) diff --git a/Detectors/TPC/workflow/include/TPCWorkflow/TPCPressureTemperatureSpec.h b/Detectors/TPC/workflow/include/TPCWorkflow/TPCPressureTemperatureSpec.h new file mode 100644 index 0000000000000..45a29e3010b33 --- /dev/null +++ b/Detectors/TPC/workflow/include/TPCWorkflow/TPCPressureTemperatureSpec.h @@ -0,0 +1,28 @@ +// 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. + +#ifndef O2_TPC_TPCPRESSURETEMPERATURE_SPEC +#define O2_TPC_TPCPRESSURETEMPERATURE_SPEC + +#include "Framework/DataProcessorSpec.h" +#include "DetectorsBase/Propagator.h" + +namespace o2 +{ +namespace tpc +{ + +o2::framework::DataProcessorSpec getTPCPressureTemperatureSpec(); + +} // end namespace tpc +} // end namespace o2 + +#endif diff --git a/Detectors/TPC/workflow/src/TPCPressureTemperatureSpec.cxx b/Detectors/TPC/workflow/src/TPCPressureTemperatureSpec.cxx new file mode 100644 index 0000000000000..e03a0ffe4308b --- /dev/null +++ b/Detectors/TPC/workflow/src/TPCPressureTemperatureSpec.cxx @@ -0,0 +1,115 @@ +// 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 TPCPressureTemperatureSpec.cxx +/// \brief device for providing pressure and temperature values +/// \author Matthias Kleiner +/// \date Jun 4, 2025 + +#include "TPCWorkflow/TPCPressureTemperatureSpec.h" +#include "Framework/Task.h" +#include "Framework/DataProcessorSpec.h" +#include "Framework/ConfigParamRegistry.h" +#include "DetectorsBase/GRPGeomHelper.h" +#include "CommonUtils/TreeStreamRedirector.h" +#include "TPCCalibration/PressureTemperatureHelper.h" + +using namespace o2::framework; + +namespace o2 +{ +namespace tpc +{ + +class PressureTemperatureDevice : public o2::framework::Task +{ + public: + PressureTemperatureDevice(std::shared_ptr req) : mCCDBRequest(req){}; + void init(o2::framework::InitContext& ic) final + { + o2::base::GRPGeomHelper::instance().setRequest(mCCDBRequest); + const int intInterval = ic.options().get("fit-interval"); + mPTHelper.setFitIntervalTemp(intInterval * 1000); + const bool enableDebugTree = ic.options().get("enable-root-output"); + if (enableDebugTree) { + mStreamer = std::make_unique("pt.root", "recreate"); + } + }; + + void endOfStream(EndOfStreamContext& eos) final + { + if (mStreamer) { + mStreamer->Close(); + } + } + + void run(o2::framework::ProcessingContext& pc) final + { + o2::base::GRPGeomHelper::instance().checkUpdates(pc); + mPTHelper.extractCCDBInputs(pc); + const auto orbitResetTimeMS = o2::base::GRPGeomHelper::instance().getOrbitResetTimeMS(); + const auto firstTFOrbit = pc.services().get().firstTForbit; + const uint64_t timestamp = orbitResetTimeMS + firstTFOrbit * o2::constants::lhc::LHCOrbitMUS * 0.001; + mPTHelper.sendPTForTS(pc, timestamp); + + if (mStreamer) { + const float pressure = mPTHelper.getPressure(timestamp); + const auto temp = mPTHelper.getTemperature(timestamp); + (*mStreamer) << "pt" + << "pressure=" << pressure + << "temperatureA=" << temp.first + << "temperatureC=" << temp.second + << "time=" << timestamp + << "\n"; + } + } + + void finaliseCCDB(o2::framework::ConcreteDataMatcher& matcher, void* obj) final + { + o2::base::GRPGeomHelper::instance().finaliseCCDB(matcher, obj); + mPTHelper.accountCCDBInputs(matcher, obj); + } + + private: + PressureTemperatureHelper mPTHelper; + std::shared_ptr mCCDBRequest; ///< info for CCDB request + std::unique_ptr mStreamer; ///< debug streamer +}; + +o2::framework::DataProcessorSpec getTPCPressureTemperatureSpec() +{ + std::vector inputs; + std::vector outputs; + o2::header::DataDescription dataDescription; + + PressureTemperatureHelper::requestCCDBInputs(inputs); + PressureTemperatureHelper::setOutputs(outputs); + + auto ccdbRequest = std::make_shared(true, // orbitResetTime + false, // GRPECS=true for nHBF per TF + false, // GRPLHCIF + false, // GRPMagField + false, // askMatLUT + o2::base::GRPGeomRequest::None, // geometry + inputs); + return DataProcessorSpec{ + "tpc-pressure-temperature", + inputs, + outputs, + AlgorithmSpec{adaptFromTask(ccdbRequest)}, + Options{ + {"enable-root-output", VariantType::Bool, false, {"Enable root-files output writers"}}, + {"fit-interval", VariantType::Int, 300, {"interval in seconds for which to e.g. perform fits of the temperature sensors"}}} // end Options + }; +} + +} // end namespace tpc +} // end namespace o2 diff --git a/Detectors/TPC/workflow/src/tpc-pressure-temperature.cxx b/Detectors/TPC/workflow/src/tpc-pressure-temperature.cxx new file mode 100644 index 0000000000000..a56067c9d843c --- /dev/null +++ b/Detectors/TPC/workflow/src/tpc-pressure-temperature.cxx @@ -0,0 +1,36 @@ +// 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 tpc-pressure-temperature.cxx +/// \author Matthias Kleiner, mkleiner@ikf.uni-frankfurt.de + +#include "TPCWorkflow/TPCPressureTemperatureSpec.h" +#include "CommonUtils/ConfigurableParam.h" +#include "Framework/ConfigParamSpec.h" + +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::swap(workflowOptions, options); +} + +#include "Framework/runDataProcessing.h" + +WorkflowSpec defineDataProcessing(ConfigContext const& config) +{ + WorkflowSpec workflow; + o2::conf::ConfigurableParam::updateFromString(config.options().get("configKeyValues")); + workflow.emplace_back(o2::tpc::getTPCPressureTemperatureSpec()); + return workflow; +}