19#include "TLinearFitter.h"
28 {
"TPC_PT_351_TEMPERATURE", 0},
29 {
"TPC_PT_376_TEMPERATURE", 1},
30 {
"TPC_PT_415_TEMPERATURE", 2},
31 {
"TPC_PT_447_TEMPERATURE", 3},
32 {
"TPC_PT_477_TEMPERATURE", 4},
33 {
"TPC_PT_488_TEMPERATURE", 5},
34 {
"TPC_PT_537_TEMPERATURE", 6},
35 {
"TPC_PT_575_TEMPERATURE", 7},
36 {
"TPC_PT_589_TEMPERATURE", 8},
37 {
"TPC_PT_629_TEMPERATURE", 9},
38 {
"TPC_PT_664_TEMPERATURE", 10},
39 {
"TPC_PT_695_TEMPERATURE", 11},
40 {
"TPC_PT_735_TEMPERATURE", 12},
41 {
"TPC_PT_757_TEMPERATURE", 13},
42 {
"TPC_PT_797_TEMPERATURE", 14},
43 {
"TPC_PT_831_TEMPERATURE", 15},
44 {
"TPC_PT_851_TEMPERATURE", 16},
45 {
"TPC_PT_895_TEMPERATURE", 17},
50 for (
size_t i = 0;
i <
raw.size(); ++
i) {
51 raw[
i].sensorNumber =
i;
66 for (
size_t i = 0;
i <
states.size(); ++
i) {
73 {StackState::OFF,
"OFF"},
74 {StackState::STBY_CONFIGURED,
"STBY_CONFIGURED"},
75 {StackState::INTERMEDIATE,
"INTERMEDIATE"},
76 {StackState::ON,
"ON"},
77 {StackState::ERROR,
"ERROR"},
78 {StackState::ERROR_LOCAL,
"ERROR_LOCAL"},
79 {StackState::SOFT_INTERLOCK,
"SOFT_INTERLOCK"},
80 {StackState::INTERLOCK,
"INTERLOCK"},
81 {StackState::RAMPIG_UP_LOW,
"RAMPIG_UP_LOW"},
82 {StackState::RAMPIG_DOWN_LOW,
"RAMPIG_DOWN_LOW"},
83 {StackState::RAMPIG_UP,
"RAMPIG_UP"},
84 {StackState::RAMPIG_DOWN,
"RAMPIG_DOWN"},
85 {StackState::MIXED,
"MIXED"},
86 {StackState::NO_CONTROL,
"NO_CONTROL"},
91 constexpr auto max = std::numeric_limits<dcs::TimeStampType>::max();
92 const std::vector<TimeStampType>
times{
103 return *std::min_element(
times.begin(),
times.end());
108 constexpr auto min = 0;
109 const std::vector<TimeStampType>
times{
120 return *std::max_element(
times.begin(),
times.end());
125 constexpr auto max = std::numeric_limits<dcs::TimeStampType>::max();
126 const std::vector<TimeStampType>
times{
132 return *std::min_element(
times.begin(),
times.end());
137 constexpr auto min = 0;
138 const std::vector<TimeStampType>
times{
144 return *std::max_element(
times.begin(),
times.end());
147bool Temperature::makeFit(TLinearFitter& fitter,
const int nDim, std::vector<double>& xVals, std::vector<double>& temperatures)
149 const int minPointsForFit = 5;
150 if (temperatures.empty() || (temperatures.size() < minPointsForFit)) {
151 LOGP(warning,
"Number of points {} for fit smaller than minimum of {}!", temperatures.size(), minPointsForFit);
155 fitter.ClearPoints();
156 fitter.AssignData(temperatures.size(), nDim, xVals.data(), temperatures.data());
157 int status = fitter.Eval();
159 LOGP(warning,
"Fit failed!");
173 TLinearFitter fitter(nDim,
"1 ++ x0 ++ x1",
"");
174 std::array<size_t, dcs::Temperature::SensorsPerSide> startPos{};
181 const int lastIntervalDuration = (refTimeMax - refTime) % fitInterval;
184 const bool procLastInt = (lastIntervalDuration / fitInterval > 0.5);
185 int numIntervals = (refTimeMax - refTime) / fitInterval + procLastInt;
186 if (numIntervals == 0) {
191 std::vector<double> xVals;
192 std::vector<double> temperatures;
193 xVals.reserve(2 * 1000);
194 temperatures.reserve(1000);
196 for (
int interval = 0; interval < numIntervals; ++interval) {
201 temperatures.clear();
208 const auto& sensor =
raw[iSensor + sensorOffset];
210 LOGP(
debug,
"sensor {}, start {}, size {}", sensor.sensorNumber, startPos[iSensor], sensor.data.size());
211 while (startPos[iSensor] < sensor.data.size()) {
212 const auto& dataPoint = sensor.data[startPos[iSensor]];
213 if (((dataPoint.time - timeStart) >= fitInterval) && (interval != numIntervals - 1)) {
214 LOGP(
debug,
"sensor {}, {} - {} >= {}", sensor.sensorNumber, dataPoint.time, timeStart, fitInterval);
217 firstTime = std::min(firstTime, dataPoint.time);
218 LastTime = std::max(LastTime, dataPoint.time);
219 const auto temperature = dataPoint.value;
222 if (temperature < 15 || temperature > 25) {
226 xVals.emplace_back(
pos.x);
227 xVals.emplace_back(
pos.y);
228 temperatures.emplace_back(temperature);
231 if (firstTime < std::numeric_limits<dcs::TimeStampType>::max() && !temperatures.empty()) {
232 const bool fitOk = makeFit(fitter, nDim, xVals, temperatures);
236 auto& stat = stats.data.emplace_back();
237 stat.time = (firstTime + LastTime) / 2;
238 stat.value.mean = fitter.GetParameter(0);
239 stat.value.gradX = fitter.GetParameter(1);
240 stat.value.gradY = fitter.GetParameter(2);
243 const float maxDeltaT = 1;
244 const float meanTemp = fitter.GetParameter(0);
245 const bool isDataGood = std::all_of(temperatures.begin(), temperatures.end(), [meanTemp, maxDeltaT](
double t) { return std::abs(t - meanTemp) < maxDeltaT; });
249 std::vector<double> xVals2;
250 std::vector<double> temperatures2;
251 xVals2.reserve(xVals.size());
252 temperatures2.reserve(temperatures.size());
253 for (
int i = 0;
i < temperatures.size(); ++
i) {
254 if (std::abs(temperatures[
i] - meanTemp) < maxDeltaT) {
255 const int idx = 2 *
i;
256 xVals2.emplace_back(xVals[idx]);
257 xVals2.emplace_back(xVals[idx + 1]);
258 temperatures2.emplace_back(temperatures[
i]);
261 const bool fitOk2 = makeFit(fitter, nDim, xVals2, temperatures2);
263 stat.value.mean = fitter.GetParameter(0);
264 stat.value.gradX = fitter.GetParameter(1);
265 stat.value.gradY = fitter.GetParameter(2);
274 if (sensor ==
"CavernAtmosPressure") {
276 }
else if (sensor ==
"CavernAtmosPressure2") {
278 }
else if (sensor ==
"SurfaceAtmosPressure") {
281 LOGP(warning,
"Unknown pressure sensor {}", sensor);
291 auto removeOutliers = [](
auto& dataVec,
auto minVal,
auto maxVal) {
293 std::remove_if(dataVec.begin(), dataVec.end(),
294 [minVal, maxVal](
const auto& dp) {
295 return (dp.value < minVal || dp.value > maxVal);
320void fillBuffer(std::pair<std::vector<float>, std::vector<TimeStampType>>&
buffer,
const std::pair<std::vector<float>, std::vector<TimeStampType>>&
values,
TimeStampType tStart,
const int minPoints)
322 const auto itStartBuff = std::lower_bound(
buffer.second.begin(),
buffer.second.end(), tStart);
323 size_t idxStartBuffer = std::distance(
buffer.second.begin(), itStartBuff);
324 if (
buffer.first.size() - idxStartBuffer < minPoints) {
325 if (
buffer.first.size() < minPoints) {
328 idxStartBuffer =
buffer.first.size() - minPoints;
332 std::pair<std::vector<float>, std::vector<TimeStampType>> buffTmp{
333 std::vector<float>(
buffer.first.begin() + idxStartBuffer,
buffer.first.end()),
334 std::vector<TimeStampType>(
buffer.second.begin() + idxStartBuffer,
buffer.second.end())};
336 buffTmp.first.insert(buffTmp.first.end(),
values.first.begin(),
values.first.end());
337 buffTmp.second.insert(buffTmp.second.end(),
values.second.begin(),
values.second.end());
339 buffer = std::move(buffTmp);
349 tStart = tStart / 1000 * 1000;
351 const int minPointsRef = 50;
356 int nIntervals = std::round((tEnd - tStart) / timeInterval);
357 if (nIntervals == 0) {
360 std::vector<TimeStampType>
times;
361 times.reserve(nIntervals);
362 for (
int i = 0;
i < nIntervals; ++
i) {
363 times.emplace_back(tStart + (
i + 0.5) * timeInterval);
367 const int minPoints = 4;
373 std::pair<std::vector<float>, std::vector<TimeStampType>> cavernAtmosPressure12;
374 std::pair<std::vector<float>, std::vector<TimeStampType>> cavernAtmosPressure1S;
375 std::pair<std::vector<float>, std::vector<TimeStampType>> cavernAtmosPressure2S;
376 cavernAtmosPressure12.first.reserve(nIntervals);
377 cavernAtmosPressure1S.first.reserve(nIntervals);
378 cavernAtmosPressure2S.first.reserve(nIntervals);
379 cavernAtmosPressure12.second.reserve(nIntervals);
380 cavernAtmosPressure1S.second.reserve(nIntervals);
381 cavernAtmosPressure2S.second.reserve(nIntervals);
383 for (
int i = 0;
i < nIntervals;
i++) {
385 const int maxDist = 600 * 1000;
386 const bool cavernOk = (cavernAtmosPressureStats.median[
i] > 0) && (cavernAtmosPressureStats.closestDistanceL[
i] < maxDist) && (cavernAtmosPressureStats.closestDistanceR[
i] < maxDist);
387 const bool cavern2Ok = (cavernAtmosPressure2Stats.median[
i] > 0) && (cavernAtmosPressure2Stats.closestDistanceL[
i] < maxDist) && (cavernAtmosPressure2Stats.closestDistanceR[
i] < maxDist);
388 const bool surfaceOk = (surfaceAtmosPressureStats.median[
i] > 0) && (surfaceAtmosPressureStats.closestDistanceL[
i] < maxDist) && (surfaceAtmosPressureStats.closestDistanceR[
i] < maxDist);
390 if (cavernOk && cavern2Ok) {
391 cavernAtmosPressure12.first.emplace_back(cavernAtmosPressureStats.median[
i] - cavernAtmosPressure2Stats.median[
i]);
392 cavernAtmosPressure12.second.emplace_back(
times[
i]);
394 if (cavernOk && surfaceOk) {
395 cavernAtmosPressure1S.first.emplace_back(cavernAtmosPressureStats.median[
i] - surfaceAtmosPressureStats.median[
i]);
396 cavernAtmosPressure1S.second.emplace_back(
times[
i]);
398 if (cavern2Ok && surfaceOk) {
399 cavernAtmosPressure2S.first.emplace_back(cavernAtmosPressure2Stats.median[
i] - surfaceAtmosPressureStats.median[
i]);
400 cavernAtmosPressure2S.second.emplace_back(
times[
i]);
414 const float maxDist = 20 * timeInterval;
415 const float maxDiff = 0.2;
416 std::pair<std::vector<float>, std::vector<TimeStampType>> robustPressureTmp;
417 robustPressureTmp.first.reserve(nIntervals);
418 robustPressureTmp.second.reserve(nIntervals);
419 std::vector<uint8_t> isOk(nIntervals);
421 for (
int i = 0;
i < nIntervals; ++
i) {
423 const float delta12 = cavernAtmosPressureStats.median[
i] - cavernAtmosPressure2Stats.median[
i] - cavernAtmosPressure12Stats.median[
i];
424 const float delta1S = cavernAtmosPressureStats.median[
i] - surfaceAtmosPressureStats.median[
i] - cavernAtmosPressure1SStats.median[
i];
425 const float delta2S = cavernAtmosPressure2Stats.median[
i] - surfaceAtmosPressureStats.median[
i] - cavernAtmosPressure2SStats.median[
i];
427 const auto distCavernAtmosPressureL = cavernAtmosPressureStats.closestDistanceL[
i];
428 const auto distCavernAtmosPressure2L = cavernAtmosPressure2Stats.closestDistanceL[
i];
429 const auto distSurfaceAtmosPressureL = surfaceAtmosPressureStats.closestDistanceL[
i];
430 const auto distCavernAtmosPressureR = cavernAtmosPressureStats.closestDistanceR[
i];
431 const auto distCavernAtmosPressure2R = cavernAtmosPressure2Stats.closestDistanceR[
i];
432 const auto distSurfaceAtmosPressureR = surfaceAtmosPressureStats.closestDistanceR[
i];
435 const bool cavernDistOk = (cavernAtmosPressureStats.median[
i] > 0) && ((distCavernAtmosPressureL < maxDist) || (distCavernAtmosPressureR < maxDist));
436 const bool cavern2DistOk = (cavernAtmosPressure2Stats.median[
i] > 0) && ((distCavernAtmosPressure2L < maxDist) || (distCavernAtmosPressure2R < maxDist));
437 const bool surfaceDistOk = (surfaceAtmosPressureStats.median[
i] > 0) && ((distSurfaceAtmosPressureL < maxDist) || (distSurfaceAtmosPressureR < maxDist));
438 const bool onlyOneSensor = (cavernDistOk + cavern2DistOk + surfaceDistOk) == 1;
440 uint8_t maskIsOkTmp = 0;
441 const int cavernBit = 0;
442 const int cavern2Bit = 1;
443 const int surfaceBit = 2;
447 if (((std::abs(delta12) < maxDiff) && (cavernDistOk && cavern2DistOk)) || onlyOneSensor) {
449 maskIsOkTmp |= (1 << cavernBit);
452 maskIsOkTmp |= (1 << cavern2Bit);
458 if ((std::abs(delta1S) < maxDiff) && ((cavernDistOk && surfaceDistOk)) || onlyOneSensor) {
460 maskIsOkTmp |= (1 << cavernBit);
463 maskIsOkTmp |= (1 << surfaceBit);
469 if ((std::abs(delta2S) < maxDiff) && ((cavern2DistOk && surfaceDistOk)) || onlyOneSensor) {
471 maskIsOkTmp |= (1 << cavern2Bit);
474 maskIsOkTmp |= (1 << surfaceBit);
480 int pressureCount = 0;
481 if ((maskIsOkTmp >> cavernBit) & 1) {
482 pressure += cavernAtmosPressureStats.median[
i];
486 if ((maskIsOkTmp >> cavern2Bit) & 1) {
487 pressure += cavernAtmosPressure2Stats.median[
i] + cavernAtmosPressure12Stats.median[
i];
491 if ((maskIsOkTmp >> surfaceBit) & 1) {
492 pressure += surfaceAtmosPressureStats.median[
i] + cavernAtmosPressure1SStats.median[
i];
496 isOk[
i] = maskIsOkTmp;
497 if (pressureCount > 0) {
498 pressure /= pressureCount;
499 robustPressureTmp.first.emplace_back(pressure);
500 robustPressureTmp.second.emplace_back(
times[
i]);
513 pOut.
isOk = std::move(isOk);
524 tree->SetAlias(
"cavernDistOk",
"robustPressure.cavernAtmosPressure.median>0 && (robustPressure.cavernAtmosPressure.closestDistanceR<robustPressure.maxDist || robustPressure.cavernAtmosPressure.closestDistanceL<robustPressure.maxDist)");
525 tree->SetAlias(
"cavern2DistOk",
"robustPressure.cavernAtmosPressure2.median>0 && (robustPressure.cavernAtmosPressure2.closestDistanceR<robustPressure.maxDist || robustPressure.cavernAtmosPressure2.closestDistanceL<robustPressure.maxDist)");
526 tree->SetAlias(
"surfaceDistOk",
"robustPressure.surfaceAtmosPressure.median>0 && (robustPressure.surfaceAtmosPressure.closestDistanceR<robustPressure.maxDist || robustPressure.surfaceAtmosPressure.closestDistanceL<robustPressure.maxDist)");
527 tree->SetAlias(
"onlyOneSensor",
"(cavernDistOk + cavern2DistOk + surfaceDistOk) == 1");
528 tree->SetAlias(
"delta12",
"robustPressure.cavernAtmosPressure.median - robustPressure.cavernAtmosPressure2.median - robustPressure.cavernAtmosPressure12.median");
529 tree->SetAlias(
"delta1S",
"robustPressure.cavernAtmosPressure.median - robustPressure.surfaceAtmosPressure.median - robustPressure.cavernAtmosPressure1S.median");
530 tree->SetAlias(
"delta2S",
"robustPressure.surfaceAtmosPressure.median - robustPressure.cavernAtmosPressure2.median - robustPressure.cavernAtmosPressure2S.median");
531 tree->SetAlias(
"delta12_Ok",
"abs(delta12)<robustPressure.maxDiff");
532 tree->SetAlias(
"delta1S_Ok",
"abs(delta1S)<robustPressure.maxDiff");
533 tree->SetAlias(
"delta2S_Ok",
"abs(delta2S)<robustPressure.maxDiff");
std::vector< unsigned long > times
void fillBuffer(std::pair< std::vector< float >, std::vector< TimeStampType > > &buffer, const std::pair< std::vector< float >, std::vector< TimeStampType > > &values, TimeStampType tStart, const int minPoints)
DCS data point data formats.
GLsizei const GLfloat * value
GLenum GLsizei GLsizei GLint * values
RollingStats getRollingStatistics(const DataTimeType &timeData, const DataType &data, const DataTime ×, const double deltaMax, const int mNthreads, const size_t minPoints=4, const size_t nClosestPoints=4)
calculates the rolling statistics of the input data
dcs::TimeStampType getMinTime(const std::vector< dcs::DataPointVector< T > > &data, const bool roundToInterval, dcs::TimeStampType fitInterval)
dcs::TimeStampType getMaxTime(const std::vector< dcs::DataPointVector< T > > &data)
constexpr unsigned char SECTORSPERSIDE
constexpr unsigned short GEMSPERSTACK
constexpr unsigned char SIDES
constexpr unsigned short GEMSTACKSPERSECTOR
a couple of static helper functions to create timestamp values for CCDB queries or override obsolete ...
std::vector< float > median
median of rolling data
std::vector< DPType > data
void fill(const TimeStampType time, const T &value)
auto getPairOfVector() const
convert data points to a vector of pairs: pair.first -> data and pair.second -> time
void append(const DataPointVector< T > &other)
RawDPsF h2o
H2O measurement from gas chromatograph.
RawDPsF neon
neon measurement from gas chromatograph
TimeStampType getMaxTime() const
TimeStampType getMinTime() const
RawDPsF co2
CO2 measurement from gas chromatograph.
RawDPsF h2oSensor
O2 measurement from dedicated gas sensor.
RawDPsF argon
argon measurement from gas chromatograph
RawDPsF o2Sensor
O2 measurement from dedicated gas sensor.
RawDPsF n2
neon measurement from gas chromatograph
std::vector< RawDPsF > voltages
voltages per GEM stack, counting is IROCs GEM1 top, bottom, GEM2 top, bottom, .. O1 ....
static const std::unordered_map< StackState, std::string > StackStateNameMap
map state to string
std::vector< RawDPsF > currents
currents per GEM stack, counting is IROCs GEM1 top, bottom, GEM2 top, bottom, .. O1 ....
std::pair< std::vector< float >, std::vector< TimeStampType > > mPressure12Buff
! buffer for normalizing the pressure cavern 1 / cavern 2
void append(const Pressure &other)
append other pressure values
void sortAndClean(float pMin=800, float pMax=1100)
RawDPsF cavernAtmosPressure2
raw pressure in the cavern from sensor 2
RobustPressure robustPressure
combined robust pressure estimator from all three sensors
std::pair< std::vector< float >, std::vector< TimeStampType > > mPressure2SBuff
! buffer for normalizing the pressure cavern 2 / surface
RawDPsF surfaceAtmosPressure
raw pressure at the surface
void fill(std::string_view sensor, const TimeStampType time, const DataType value)
fill pressure data
static void setAliases(TTree *tree)
set aliases for the cuts used in the calculation of the robust pressure
void makeRobustPressure(TimeStampType timeInterval=100 *1000, TimeStampType timeIntervalRef=24 *60 *1000, TimeStampType tStart=1, TimeStampType tEnd=0, const int nthreads=1)
average pressure values for given time interval
TimeStampType getMinTime() const
std::pair< std::vector< float >, std::vector< TimeStampType > > mCavernAtmosPressure2Buff
! buffer for the pressure cavern 2 sensor
std::pair< std::vector< float >, std::vector< TimeStampType > > mSurfaceAtmosPressureBuff
! buffer for the pressure surface sensort
std::pair< std::vector< float >, std::vector< TimeStampType > > mPressure1SBuff
! buffer for normalizing the pressure cavern 1 / surface
RawDPsF cavernAtmosPressure
raw pressure in the cavern from sensor 1
std::pair< std::vector< float >, std::vector< TimeStampType > > mRobPressureBuff
! buffer for the robust pressure
std::pair< std::vector< float >, std::vector< TimeStampType > > mCavernAtmosPressure1Buff
! buffer for the pressure cavern 1 sensor
void clear()
\clear all stored data except the buffer
TimeStampType getMaxTime() const
float maxDist
maximum allowed time distance between sensors to be accepted for robust pressure calculation
Stats cavernAtmosPressure2S
rolling statistics of cavernAtmosPressure2/surfaceAtmosPressure
Stats cavernAtmosPressure1S
rolling statistics of cavernAtmosPressure/surfaceAtmosPressure
std::vector< uint8_t > isOk
bit mask of valid sensors: cavernBit 0, cavern2Bit = 1, surfaceBit = 2
Stats cavernAtmosPressure2
rolling statistics of cavern sensor 2
std::vector< float > robustPressure
combined robust pressure value that should be used
Stats cavernAtmosPressure12
rolling statistics of cavernAtmosPressure/cavernAtmosPressure2
float maxDiff
maximum allowed pressure difference between sensors to be accepted for robust pressure calculation
std::vector< TimeStampType > time
time stamps of all pressure values
Stats cavernAtmosPressure
rolling statistics of cavern sensor 1
TimeStampType timeIntervalRef
reference time interval used for normalization of pressure sensors
Stats surfaceAtmosPressure
rolling statistics of surface sensor
TimeStampType timeInterval
time interval used for rolling statistics
static constexpr std::array< Position, SensorsPerSide *SIDES > SensorPosition
StatsDPs statsA
statistics fit values per integration time A-Side
StatsDPs statsC
statistics fit values per integration time C-Side
void fitTemperature(Side side, dcs::TimeStampType fitInterval=5 *60 *1000, const bool roundToInterval=false)
make fit of the mean temperature and gradients in time intervals
static constexpr int SensorsPerSide
number of temperature sensors in the active volume per side
static const std::unordered_map< std::string, int > SensorNameMap
std::vector< RawDPsF > raw
raw temperature values from DCS for
VectorOfTObjectPtrs other
std::unique_ptr< TTree > tree((TTree *) flIn.Get(std::string(o2::base::NameConf::CTFTREENAME).c_str()))