TPC Update CheckOfSlices with Masking (#2245)

* Save WIP

* TPC Update CheckOfSlices with Masking

* Adding changes

---------

Co-authored-by: mlesch <ga94pez@mytum.de>

Author: mlesch

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<int> mMaskedPoints;
 
   double mMean = 0;
   double mStdev;

Modules/TPC/include/TPC/Utility.h

@@ -73,7 +73,7 @@ void getTimestamp(const std::string& metaInfo, std::vector<long>& timeStamps);
 /// \param limit Most recent timestamp to be processed
 std::vector<long> 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<long> 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<int>&, 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<int>& maskPoints);
+
+/// \brief Calculates mean and stddev from yValues of a TGraph. Overloaded function, actual calculation in retrieveStatistics
+/// \param values std::vector<double>& vector that contains the data points
+/// \param errors std::vector<double>& 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<double>& values, std::vector<double>& errors, bool useErrors, double& mean, double& stddevOfMean);
 } // namespace o2::quality_control_modules::tpc
 #endif // QUALITYCONTROL_TPCUTILITY_H

Modules/TPC/src/CheckOfSlices.cxx

@@ -21,6 +21,7 @@
 #include <fmt/format.h>
 #include "Common/Utils.h"
 #include "TPC/Utility.h"
+#include "Algorithm/RangeTokenizer.h"
 
 #include <TCanvas.h>
 #include <TGraphErrors.h>
@@ -32,6 +33,7 @@
 
 #include <vector>
 #include <numeric>
+#include <sstream>
 
 namespace o2::quality_control_modules::tpc
 {
@@ -98,6 +100,12 @@ void CheckOfSlices::configure()
     }
     mExpectedPhysicsValue = common::getFromConfig<double>(mCustomParameters, "expectedPhysicsValue", 1);
   }
+
+  std::string maskingValues = common::getFromConfig<string>(mCustomParameters, "maskedPoints", "");
+  mMaskedPoints.clear();
+  if (maskingValues != "") {
+    mMaskedPoints = RangeTokenizer::tokenize<int>(maskingValues);
+  }
 }
 
 Quality CheckOfSlices::check(std::map<std::string, std::shared_ptr<MonitorObject>>* moMap)
@@ -168,7 +176,7 @@ Quality CheckOfSlices::check(std::map<std::string, std::shared_ptr<MonitorObject
     }
   }
 
-  calculateStatistics(yValues, yErrors, useErrors, 0, NBins, mMean, mStdev);
+  calculateStatistics(yValues, yErrors, useErrors, 0, NBins, mMean, mStdev, mMaskedPoints);
 
   for (size_t i = 0; i < v.size(); ++i) {
 

Modules/TPC/src/Utility.cxx

@@ -30,6 +30,7 @@
 // external includes
 #include <Framework/Logger.h>
 #include <bitset>
+#include <algorithm>
 
 namespace o2::quality_control_modules::tpc
 {
@@ -210,44 +211,101 @@ 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<double> v(yValues + firstPoint, yValues + lastPoint);
+  std::vector<double> vErr;
+
+  if (useErrors) {
+    const std::vector<double> 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<int>& 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<double> v;
+  const std::vector<double> 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<double> vErr;
+  if (useErrors) {
+    const std::vector<double> 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<double>& values, std::vector<double>& 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<double> 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<double> 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;
     }
   }
 
   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 {
       stddevOfMean = sqrt(1. / sumOfWeights);
     }
   } else { // for >= 2 points, we calculate the spread
     if (!useErrors) {
-      std::vector<double> diff(v.size());
-      std::transform(v.begin(), v.end(), diff.begin(), [mean](double x) { return x - mean; });
+      std::vector<double> 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);