diff --git a/Modules/TPC/include/TPC/CheckOfSlices.h b/Modules/TPC/include/TPC/CheckOfSlices.h index e19ce1bb88..ab16896c08 100644 --- a/Modules/TPC/include/TPC/CheckOfSlices.h +++ b/Modules/TPC/include/TPC/CheckOfSlices.h @@ -55,6 +55,7 @@ class CheckOfSlices : public o2::quality_control::checker::CheckInterface double mRangeMedium; double mRangeBad; bool mSliceTrend; + std::vector mMaskedPoints; double mMean = 0; double mStdev; diff --git a/Modules/TPC/include/TPC/Utility.h b/Modules/TPC/include/TPC/Utility.h index 3478a6268d..18f505200f 100644 --- a/Modules/TPC/include/TPC/Utility.h +++ b/Modules/TPC/include/TPC/Utility.h @@ -73,7 +73,7 @@ void getTimestamp(const std::string& metaInfo, std::vector& timeStamps); /// \param limit Most recent timestamp to be processed std::vector getDataTimestamps(const o2::ccdb::CcdbApi& cdbApi, const std::string_view path, const unsigned int nFiles, const long limit); -/// \brief Calculates mean and stddev from yValues of a TGraph +/// \brief Calculates mean and stddev from yValues of a TGraph. Overloaded function, actual calculation in retrieveStatistics /// \param yValues const double* pointer to yValues of TGraph (via TGraph->GetY()) /// \param yErrors const double* pointer to y uncertainties of TGraph (via TGraph->GetEY()) /// \param useErrors bool whether uncertainties should be used in calculation of mean and stddev of mean @@ -82,5 +82,24 @@ std::vector getDataTimestamps(const o2::ccdb::CcdbApi& cdbApi, const std:: /// \param mean double&, reference to double that should store mean /// \param stddevOfMean double&, reference to double that should store stddev of mean void calculateStatistics(const double* yValues, const double* yErrors, bool useErrors, const int firstPoint, const int lastPoint, double& mean, double& stddevOfMean); + +/// \brief Calculates mean and stddev from yValues of a TGraph. Overloaded function +/// \param yValues const double* pointer to yValues of TGraph (via TGraph->GetY()) +/// \param yErrors const double* pointer to y uncertainties of TGraph (via TGraph->GetEY()) +/// \param useErrors bool whether uncertainties should be used in calculation of mean and stddev of mean +/// \param firstPoint const int, first point of yValues to include in calculation +/// \param lastPoint const int, last point of yValues to include in calculation +/// \param mean double&, reference to double that should store mean +/// \param stddevOfMean double&, reference to double that should store stddev of mean +/// \param maskPoints std::vector&, points of the selected TGraph-points that should be masked +void calculateStatistics(const double* yValues, const double* yErrors, bool useErrors, const int firstPoint, const int lastPoint, double& mean, double& stddevOfMean, std::vector& maskPoints); + +/// \brief Calculates mean and stddev from yValues of a TGraph. Overloaded function, actual calculation in retrieveStatistics +/// \param values std::vector& vector that contains the data points +/// \param errors std::vector& vector that contains the data errors +/// \param useErrors bool whether uncertainties should be used in calculation of mean and stddev of mean +/// \param mean double&, reference to double that should store mean +/// \param stddevOfMean double&, reference to double that should store stddev of mean +void retrieveStatistics(std::vector& values, std::vector& errors, bool useErrors, double& mean, double& stddevOfMean); } // namespace o2::quality_control_modules::tpc #endif // QUALITYCONTROL_TPCUTILITY_H \ No newline at end of file diff --git a/Modules/TPC/src/CheckOfSlices.cxx b/Modules/TPC/src/CheckOfSlices.cxx index b075982751..f281d39596 100644 --- a/Modules/TPC/src/CheckOfSlices.cxx +++ b/Modules/TPC/src/CheckOfSlices.cxx @@ -21,6 +21,7 @@ #include #include "Common/Utils.h" #include "TPC/Utility.h" +#include "Algorithm/RangeTokenizer.h" #include #include @@ -32,6 +33,7 @@ #include #include +#include namespace o2::quality_control_modules::tpc { @@ -98,6 +100,12 @@ void CheckOfSlices::configure() } mExpectedPhysicsValue = common::getFromConfig(mCustomParameters, "expectedPhysicsValue", 1); } + + std::string maskingValues = common::getFromConfig(mCustomParameters, "maskedPoints", ""); + mMaskedPoints.clear(); + if (maskingValues != "") { + mMaskedPoints = RangeTokenizer::tokenize(maskingValues); + } } Quality CheckOfSlices::check(std::map>* moMap) @@ -168,7 +176,7 @@ Quality CheckOfSlices::check(std::map #include +#include namespace o2::quality_control_modules::tpc { @@ -210,25 +211,82 @@ void calculateStatistics(const double* yValues, const double* yErrors, bool useE return; } + if (useErrors && !yErrors) { + ILOG(Error, Support) << "In calculateStatistics(): requested to use errors of data but TGraph does not contain errors." << ENDM; + useErrors = false; + } + + std::vector v(yValues + firstPoint, yValues + lastPoint); + std::vector vErr; + + if (useErrors) { + const std::vector vErr_temp(yErrors + firstPoint, yErrors + lastPoint); + for (int i = 0; i < vErr_temp.size(); i++) { + vErr.push_back(vErr_temp[i]); + } + } + + retrieveStatistics(v, vErr, useErrors, mean, stddevOfMean); +} + +void calculateStatistics(const double* yValues, const double* yErrors, bool useErrors, const int firstPoint, const int lastPoint, double& mean, double& stddevOfMean, std::vector& maskPoints) +{ + // yErrors returns nullptr for TGraph (no errors) + if (lastPoint - firstPoint <= 0) { + ILOG(Error, Support) << "In calculateStatistics(), the first and last point of the range have to differ!" << ENDM; + return; + } + + if (useErrors && !yErrors) { + ILOG(Error, Support) << "In calculateStatistics(): requested to use errors of data but TGraph does not contain errors." << ENDM; + useErrors = false; + } + + std::vector v; + const std::vector v_temp(yValues + firstPoint, yValues + lastPoint); + for (int i = 0; i < v_temp.size(); i++) { + if (std::find(maskPoints.begin(), maskPoints.end(), i) == maskPoints.end()) { // i is not in the masked points + v.push_back(v_temp[i]); + } + } + + std::vector vErr; + if (useErrors) { + const std::vector vErr_temp(yErrors + firstPoint, yErrors + lastPoint); + for (int i = 0; i < vErr_temp.size(); i++) { + if (std::find(maskPoints.begin(), maskPoints.end(), i) == maskPoints.end()) { // i is not in the masked points + vErr.push_back(vErr_temp[i]); + } + } + } + + retrieveStatistics(v, vErr, useErrors, mean, stddevOfMean); +} + +void retrieveStatistics(std::vector& values, std::vector& errors, bool useErrors, double& mean, double& stddevOfMean) +{ + if ((errors.size() != values.size()) && useErrors) { + ILOG(Error, Support) << "In retrieveStatistics(): errors do not match data points, omitting errors" << ENDM; + useErrors = false; + } + double sum = 0.; double sumSquare = 0.; double sumOfWeights = 0.; // sum w_i double sumOfSquaredWeights = 0.; // sum (w_i)^2 double weight = 0.; - const std::vector v(yValues + firstPoint, yValues + lastPoint); if (!useErrors) { // In case of no errors, we set our weights equal to 1 - sum = std::accumulate(v.begin(), v.end(), 0.0); - sumOfWeights = v.size(); - sumOfSquaredWeights = v.size(); + sum = std::accumulate(values.begin(), values.end(), 0.0); + sumOfWeights = values.size(); + sumOfSquaredWeights = values.size(); } else { // In case of errors, we set our weights equal to 1/sigma_i^2 - const std::vector vErr(yErrors + firstPoint, yErrors + lastPoint); - for (size_t i = 0; i < v.size(); i++) { - weight = 1. / std::pow(vErr[i], 2.); - sum += v[i] * weight; - sumSquare += v[i] * v[i] * weight; + for (size_t i = 0; i < values.size(); i++) { + weight = 1. / std::pow(errors[i], 2.); + sum += values[i] * weight; + sumSquare += values[i] * values[i] * weight; sumOfWeights += weight; sumOfSquaredWeights += weight * weight; } @@ -236,7 +294,7 @@ void calculateStatistics(const double* yValues, const double* yErrors, bool useE mean = sum / sumOfWeights; - if (v.size() == 1) { // we only have one point, we keep it's uncertainty + if (values.size() == 1) { // we only have one point, we keep it's uncertainty if (!useErrors) { stddevOfMean = 0.; } else { @@ -244,10 +302,10 @@ void calculateStatistics(const double* yValues, const double* yErrors, bool useE } } else { // for >= 2 points, we calculate the spread if (!useErrors) { - std::vector diff(v.size()); - std::transform(v.begin(), v.end(), diff.begin(), [mean](double x) { return x - mean; }); + std::vector diff(values.size()); + std::transform(values.begin(), values.end(), diff.begin(), [mean](double x) { return x - mean; }); double sq_sum = std::inner_product(diff.begin(), diff.end(), diff.begin(), 0.0); - stddevOfMean = std::sqrt(sq_sum / (v.size() * (v.size() - 1.))); + stddevOfMean = std::sqrt(sq_sum / (values.size() * (values.size() - 1.))); } else { double ratioSumWeight = sumOfSquaredWeights / (sumOfWeights * sumOfWeights); stddevOfMean = sqrt((sumSquare / sumOfWeights - mean * mean) * (1. / (1. - ratioSumWeight)) * ratioSumWeight);