general_utility.h

Go to the documentation of this file.

//  --- --- ---
//  General utility functions
//
//  Auth: Nicola Rubini
//  Mail: nicola.rubini@bo.infn.it
//
#pragma once
 
#include "graphutils.C"
 
namespace utility
{
  //  Swap two values
  template <typename arg_type>
  inline void swap_values(arg_type &first_element, arg_type &second_element);
  //  Square sum
  template <typename arg_type>
  inline arg_type sq_sum(arg_type value);
  template <typename arg_type, typename... Args>
  inline arg_type sq_sum(arg_type first, Args... args);
  //  Average
  // template <typename T, typename = typename std::enable_if<std::is_arithmetic<T>::value>::type>
  // std::map<std::string, std::array<T, 2>> average(std::vector<std::array<T, 2>> list_of_measurements, bool skip_unfit_skim = false);
  //  Set precision of the number
  template <typename T>
  T round_digits(T value, int digits);
  //  Concatenate vectors
  template <typename T>
  std::vector<T> merge(const std::vector<T> &vec1, const std::vector<T> &vec2);
  //  Get folder list of a location
  std::vector<std::string> get_folders(const std::string &dir_path);
 
  //  Read a txt file with some numbers
  void readFileToMap(const std::string &filename, std::map<std::string, std::vector<std::string>> &dataMap);
  void readTxtToMap(const std::string &filename, std::map<std::string, std::vector<std::string>> &dataMap);
  void readCsvToMap(const std::string &filename, std::map<std::string, std::vector<std::string>> &dataMap);
 
  //  Calculate intercept between lines
  std::array<std::array<float, 2>, 2> get_intercept(float x0_1, float ex0_1, float x1_1, float ex1_1, float x0_2, float ex0_2, float x1_2, float ex1_2);
  std::array<std::array<float, 2>, 2> get_intercept_pol1_pol1(float x0_1, float ex0_1, float x1_1, float ex1_1, float x0_2, float ex0_2, float x1_2, float ex1_2);
  std::array<std::array<float, 2>, 2> get_intercept(std::array<float, 2> m1, std::array<float, 2> q1, std::array<float, 2> m2, std::array<float, 2> q2) { return get_intercept(m1[0], m1[1], q1[0], q1[1], m2[0], m2[1], q2[0], q2[1]); };
  std::array<std::array<float, 2>, 2> get_intercept_x(float m, float em, float q, float eq) { return utility::get_intercept(m, em, q, eq, 0, 0, 0, 0); };
  std::array<std::array<float, 2>, 2> get_intercept_y(float m, float em, float q, float eq) { return {q, eq}; };
  std::array<std::array<float, 2>, 2> get_intercept_pol1_pol2(float x0_1, float ex0_1, float x1_1, float ex1_1, float x0_2, float ex0_2, float x1_2, float ex1_2, float x2_2, float ex2_2);
  std::array<std::array<float, 2>, 2> get_intercept_tf1(TF1 *first_target, TF2 *second_target, double first_guess, double starting_coarse, double req_precision);
}
 
template <typename arg_type>
inline void utility::swap_values(arg_type &first_element, arg_type &second_element)
{
  auto memory = second_element;
  second_element = first_element;
  first_element = memory;
}
 
template <typename arg_type>
inline arg_type utility::sq_sum(arg_type value)
{
  return value * value; // Base case: square the single argument
}
 
template <typename arg_type, typename... Args>
inline arg_type utility::sq_sum(arg_type first, Args... args)
{
  return first * first + sq_sum(args...); // Recursive call
}
 
template <typename T>
T utility::round_digits(T value, int digits)
{
  //  Check the passed argument is a number
  static_assert(std::is_arithmetic<T>::value, "Type must be numeric");
  T factor = std::pow(10, digits);
  return std::ceil(value * factor) / factor;
}
 
template <typename T>
std::vector<T> utility::merge(const std::vector<T> &vec1, const std::vector<T> &vec2)
{
  std::vector<T> result = vec1;                          // Start with the first vector
  result.insert(result.end(), vec2.begin(), vec2.end()); // Concatenate the second vector
  return result;
}
 
std::vector<std::string> utility::get_folders(const std::string &dir_path)
{
  std::vector<std::string> folders;
  for (const auto &entry : std::filesystem::directory_iterator(dir_path))
    if (entry.is_directory())
      folders.push_back(entry.path().filename().string());
  return folders;
}
 
namespace utility
{
  //  Check again
  template <typename T, typename = typename std::enable_if<std::is_arithmetic<T>::value>::type>
  std::map<std::string, std::array<T, 2>> average(std::vector<std::array<T, 2>> list_of_measurements, bool skip_unfit_skim = false)
  {
    //  Final result
    std::map<std::string, std::array<T, 2>> result;
    std::vector<std::array<T, 2>> skimmed_data;
 
    //  Skim dataset from dangerous values
    skimmed_data = list_of_measurements;
    if (!skip_unfit_skim)
    {
      skimmed_data.clear();
      for (auto current_measurement : list_of_measurements)
      {
        if (std::isnan(current_measurement[0]) || std::isnan(current_measurement[1]))
          continue;
        if (std::isinf(current_measurement[0]) || std::isinf(current_measurement[1]))
          continue;
        skimmed_data.push_back(current_measurement);
      }
    }
 
    //  Measure average and error on the average
    for (auto current_measurement : skimmed_data)
    {
      result["ave_err"][0] += current_measurement[0];
      result["ave_rms"][0] += current_measurement[0];
      result["ave_sqe"][0] += current_measurement[0];
      result["sqa_err"][0] += current_measurement[0] * current_measurement[0];
      result["sqa_rms"][0] += current_measurement[0] * current_measurement[0];
      result["sqa_sqe"][0] += current_measurement[0] * current_measurement[0];
      result["ave_err"][1] += current_measurement[1];
      // result["ave_rms"][1] += current_measurement[1];
      result["ave_sqe"][1] += current_measurement[1] * current_measurement[1];
      result["sqa_err"][1] += current_measurement[1];
      // result["sqa_rms"][1] += current_measurement[1];
      result["sqa_sqe"][1] += current_measurement[1] * current_measurement[1];
    }
    result["ave_err"][0] /= skimmed_data.size();
    result["ave_rms"][0] /= skimmed_data.size();
    result["ave_sqe"][0] /= skimmed_data.size();
    result["sqa_err"][0] /= skimmed_data.size();
    result["sqa_rms"][0] /= skimmed_data.size();
    result["sqa_sqe"][0] /= skimmed_data.size();
    result["sqa_err"][0] = sqrt(result["sqa_err"][0]);
    result["sqa_rms"][0] = sqrt(result["sqa_rms"][0]);
    result["sqa_sqe"][0] = sqrt(result["sqa_sqe"][0]);
    result["ave_err"][1] /= skimmed_data.size();
    // result["ave_rms"][1] += current_measurement[1];
    result["ave_sqe"][1] = sqrt(result["ave_sqe"][1]);
    result["sqa_err"][1] /= skimmed_data.size();
    // result["sqa_rms"][1] += current_measurement[1];
    result["sqa_sqe"][1] = sqrt(result["sqa_sqe"][1]);
 
    //  Measure RMS
 
    return result;
  }
}
 
void utility::readCsvToMap(const std::string &filename, std::map<std::string, std::vector<std::string>> &dataMap)
{
  if (!std::filesystem::exists(filename))
  {
    throw std::runtime_error("[ERROR][readFileToMap] File does not exist: " + filename);
  }
 
  std::ifstream file(filename);
  if (!file)
  {
    throw std::runtime_error("[ERROR][readFileToMap] Could not open file " + filename);
  }
 
  // Read header line
  std::string headerLine;
  if (!std::getline(file, headerLine))
  {
    throw std::runtime_error("[ERROR][readFileToMap] File " + filename + " is empty or missing header row.");
  }
 
  // Ensure the string has at least 3 characters to check
  if (filename.size() < 3)
  {
    std::cerr << "[ERROR][utility::readFileToMap] Filename too short to determine file type: " << filename << std::endl;
    return;
  }
 
  // Parse column titles
  std::vector<std::string> columnTitles;
  std::istringstream headerStream(headerLine);
  std::string title;
 
  // Split header line by commas
  size_t start = 0, end = 0;
  while ((end = headerLine.find(',', start)) != std::string::npos)
  {
    columnTitles.push_back(std::regex_replace(headerLine.substr(start, end - start), std::regex("^\\s+|\\s+$"), ""));
    start = end + 1;
  }
  columnTitles.push_back(std::regex_replace(headerLine.substr(start), std::regex("^\\s+|\\s+$"), ""));
 
  for (const auto &col : columnTitles)
  {
    if (dataMap.find(col) == dataMap.end())
    {
      dataMap[col] = {}; // Initialize vector if not already present
    }
  }
 
  // Define the allowed characters regex pattern
  const std::regex validPattern(R"(^\s*[+-]?[0-9]+(\.[0-9]*)?([eE][+-]?[0-9]+)?\s*$)");
 
  // Read remaining lines
  std::string line;
  while (std::getline(file, line))
  {
    std::vector<std::string> rowValues;
    start = 0;
 
    // Split the line by commas
    while ((end = line.find(',', start)) != std::string::npos)
    {
      rowValues.push_back(line.substr(start, end - start));
      start = end + 1;
    }
    rowValues.push_back(line.substr(start));
 
    if (rowValues.size() != columnTitles.size())
    {
      throw std::runtime_error("[ERROR][readFileToMap] Mismatched data and header columns in file " + filename);
    }
 
    bool isValidLine = true;
 
    // Check each value for allowed characters
    for (const auto &value : rowValues)
    {
      if (!std::regex_match(value, validPattern))
      {
        isValidLine = false;
        break;
      }
    }
 
    // If the line is valid, add its values to the respective columns
    if (isValidLine)
    {
      for (size_t i = 0; i < columnTitles.size(); ++i)
      {
        dataMap[columnTitles[i]].push_back(rowValues[i]);
      }
    }
    else
    {
      std::cout << "[WARNING] Skipping line for character incompatibility: " << line << std::endl;
    }
  }
 
  file.close();
}
 
void utility::readTxtToMap(const std::string &filename, std::map<std::string, std::vector<std::string>> &dataMap)
{
  if (!std::filesystem::exists(filename))
  {
    throw std::runtime_error("[ERROR][readFileToMap] File does not exist: " + filename);
  }
 
  std::ifstream file(filename);
  if (!file)
  {
    throw std::runtime_error("[ERROR][readFileToMap] Could not open file " + filename);
  }
 
  // Read header line
  std::string headerLine;
  if (!std::getline(file, headerLine))
  {
    throw std::runtime_error("[ERROR][readFileToMap] File " + filename + " is empty or missing header row.");
  }
 
  // Parse column titles
  std::istringstream headerStream(headerLine);
  std::vector<std::string> columnTitles;
  std::string title;
  while (headerStream >> title)
  {
    columnTitles.push_back(title);
    if (dataMap.find(title) == dataMap.end())
    {
      dataMap[title] = {}; // Initialize vector if not already present
    }
  }
 
  // Define the allowed characters regex pattern
  const std::regex validPattern(R"(^\s*[+-]?[0-9]+(\.[0-9]*)?([eE][+-]?[0-9]+)?\s*$)");
 
  // Read remaining lines
  std::string line;
  while (std::getline(file, line))
  {
    std::istringstream lineStream(line);
    std::vector<std::string> rowValues;
    std::string value;
    bool isValidLine = true;
 
    // Process each column in the line
    for (size_t i = 0; i < columnTitles.size(); ++i)
    {
      if (!(lineStream >> value))
      {
        throw std::runtime_error("[ERROR][readFileToMap] Mismatched data and header columns.");
      }
 
      // Check if the value contains only allowed characters
      if (!std::regex_match(value, validPattern))
      {
        isValidLine = false;
        break; // Skip the rest of the line
      }
 
      rowValues.push_back(value);
    }
 
    // If the line is valid, add its values to the respective columns
    if (isValidLine)
    {
      for (size_t i = 0; i < rowValues.size(); ++i)
      {
        dataMap[columnTitles[i]].push_back(rowValues[i]);
      }
    }
    else
      cout << "[WARNING] Skipping line for character incompatibility" << endl;
  }
 
  file.close();
}
 
void utility::readFileToMap(const std::string &filename, std::map<std::string, std::vector<std::string>> &dataMap)
{
  if (!std::filesystem::exists(filename))
  {
    throw std::runtime_error("[ERROR][readFileToMap] File does not exist: " + filename);
  }
 
  std::ifstream file(filename);
  if (!file)
  {
    throw std::runtime_error("[ERROR][readFileToMap] Could not open file " + filename);
  }
 
  // Read header line
  std::string headerLine;
  if (!std::getline(file, headerLine))
  {
    throw std::runtime_error("[ERROR][readFileToMap] File " + filename + " is empty or missing header row.");
  }
 
  // Ensure the string has at least 3 characters to check
  if (filename.size() < 3)
  {
    std::cerr << "[ERROR][utility::readFileToMap] Filename too short to determine file type: " << filename << std::endl;
    return;
  }
 
  // Get the last three characters of the string
  std::string fileExtension = filename.substr(filename.size() - 3);
 
  // Convert to lowercase for case-insensitive comparison
  for (char &c : fileExtension)
    c = std::tolower(c);
 
  // Perform actions based on the file extension
  if (fileExtension == "txt")
  {
    utility::readTxtToMap(filename, dataMap);
  }
  else if (fileExtension == "csv")
  {
    utility::readCsvToMap(filename, dataMap);
  }
  else
  {
    std::cerr << "[WARNING][utility::readFileToMap] Unknown file type: " << filename << ", treating as txt" << std::endl;
  }
}
 
std::array<std::array<float, 2>, 2> utility::get_intercept(float x0_1, float ex0_1, float x1_1, float ex1_1, float x0_2, float ex0_2, float x1_2, float ex1_2)
{
  //  Intercept of pol1 and pol1:
  //  The analytic solution with the subtraction of the linear equation to the second linear equation
 
  //  Define the result container
  std::array<std::array<float, 2>, 2> result;
  return result;
}
 
std::array<std::array<float, 2>, 2> utility::get_intercept_pol1_pol1(float x0_1, float ex0_1, float x1_1, float ex1_1, float x0_2, float ex0_2, float x1_2, float ex1_2)
{
  //  Intercept of pol1 and pol1:
  //  The analytic solution with the subtraction of the linear equation to the second linear equation
 
  //  Define the result container
  std::array<std::array<float, 2>, 2> result = {0};
  return result;
 
  /*
    //  Analytic solution
    //  xi_j: (i) is the coefficient of x^i and (j) means it is related to the equation of position 1,2
    result[0][0] = (x0_2 - x0_1) / (x1_1 - x1_2);
    result[1][0] = (x1_1 * x0_2 - x1_2 * x0_1) / (x1_1 - x1_2);
 
    //  Error evaluation
    //  eabs_i_xi_j : (i) is the x (0) or y (1) relative to the xi_j component
    //  eabs: absolute error
    //  Errors are calculated with the derivative method
    auto eabs_1_x0_1 = ex0_1 * fabs(1. / (d_delta));
    auto eabs_0_x0_1 = ex0_1 * fabs(1 + x1_1 / (d_delta));
    auto eabs_1_x1_1 = ex1_1 * fabs(1 - (x1_2 - x1_1) / (2 * x2_2 * d_delta));
    auto eabs_0_x1_1 = ex1_1 * fabs((-(x1_2 - x1_1) * x1_1 / (2 * x2_2 * d_delta)) + ((d_delta) / (2 * x2_2)) + 2 * x1_1 - x1_2);
    auto eabs_1_x0_2 = ex0_2 * fabs(-1. / (d_delta));
    auto eabs_0_x0_2 = ex0_2 * fabs(-x1_1 / (d_delta));
    auto eabs_1_x1_2 = ex1_2 * fabs((x1_2 - x1_1) / (d_delta * 2 * x2_2) - 1);
    auto eabs_0_x1_2 = ex1_2 * fabs(x1_1 * (x1_2 - x1_1) / (d_delta * 2 * x2_2) - 1);
 
    //  Assign errors on result
    result[0][1] = sqrt(eabs_0_x0_1 * eabs_0_x0_1 + eabs_0_x1_1 * eabs_0_x1_1 + eabs_0_x0_2 * eabs_0_x0_2 + eabs_0_x1_2 * eabs_0_x1_2 + eabs_0_x2_2 * eabs_0_x2_2);
    result[1][1] = sqrt(eabs_1_x0_1 * eabs_1_x0_1 + eabs_1_x1_1 * eabs_1_x1_1 + eabs_1_x0_2 * eabs_1_x0_2 + eabs_1_x1_2 * eabs_1_x1_2 + eabs_1_x2_2 * eabs_1_x2_2);
    return result;
  }
 
  std::array<std::array<float, 2>, 2> utility::get_intercept_pol1_pol2(float x0_1, float ex0_1, float x1_1, float ex1_1, float x0_2, float ex0_2, float x1_2, float ex1_2, float x2_2, float ex2_2)
  {
    //  Intercept of pol1 and pol2:
    //  The analytic solution with the subtraction of the linear equation to the quadratic equation
 
    //  Define the result container
    std::array<std::array<float, 2>, 2> result;
 
    //  Analytic solution
    //  xi_j: (i) is the coefficient of x^i and (j) means it is related to the equation of degree 1,2
    result[0][0] = (-(x1_2 - x1_1) + sqrt((x1_2 - x1_1) * (x1_2 - x1_1) - 4 * x2_2 * (x0_2 - x0_1))) / (2 * x2_2);
    result[1][0] = x0_1 + result[0][0] * x1_1;
 
    //  Error evaluation
    //  eabs_i_xi_j : (i) is the x (0) or y (1) relative to the xi_j component
    //  eabs: absolute error
    //  Errors are calculated with the derivative method
    auto d_delta = sqrt((x1_2 - x1_1) * (x1_2 - x1_1) + 4 * x2_2 * (x0_1 - x0_2));
    auto eabs_1_x0_1 = ex0_1 * fabs(1. / (d_delta));
    auto eabs_0_x0_1 = ex0_1 * fabs(1 + x1_1 / (d_delta));
    auto eabs_1_x1_1 = ex1_1 * fabs(1 - (x1_2 - x1_1) / (2 * x2_2 * d_delta));
    auto eabs_0_x1_1 = ex1_1 * fabs((-(x1_2 - x1_1) * x1_1 / (2 * x2_2 * d_delta)) + ((d_delta) / (2 * x2_2)) + 2 * x1_1 - x1_2);
    auto eabs_1_x0_2 = ex0_2 * fabs(-1. / (d_delta));
    auto eabs_0_x0_2 = ex0_2 * fabs(-x1_1 / (d_delta));
    auto eabs_1_x1_2 = ex1_2 * fabs((x1_2 - x1_1) / (d_delta * 2 * x2_2) - 1);
    auto eabs_0_x1_2 = ex1_2 * fabs(x1_1 * (x1_2 - x1_1) / (d_delta * 2 * x2_2) - 1);
    auto eabs_1_x2_2 = ex2_2 * fabs(((x1_2 - x1_1) * (x1_2 - x1_1) + 2 * x2_2 * (x0_1 - x0_2)) / (d_delta * 2 * x2_2) - 1);
    auto eabs_0_x2_2 = ex2_2 * fabs(-x1_1 * ((x1_2 - x1_1) * (x1_2 - x1_1) + 2 * x2_2 * (x0_1 - x0_2)) / (d_delta * 2 * x2_2) - 1);
 
    //  Assign errors on result
    result[0][1] = sqrt(eabs_0_x0_1 * eabs_0_x0_1 + eabs_0_x1_1 * eabs_0_x1_1 + eabs_0_x0_2 * eabs_0_x0_2 + eabs_0_x1_2 * eabs_0_x1_2 + eabs_0_x2_2 * eabs_0_x2_2);
    result[1][1] = sqrt(eabs_1_x0_1 * eabs_1_x0_1 + eabs_1_x1_1 * eabs_1_x1_1 + eabs_1_x0_2 * eabs_1_x0_2 + eabs_1_x1_2 * eabs_1_x1_2 + eabs_1_x2_2 * eabs_1_x2_2);
    return result;
    */
}
 
std::array<std::array<float, 2>, 2> utility::get_intercept_tf1(TF1 *first_target, TF2 *second_target, double first_guess, double starting_coarse, double req_precision)
{
  //  Result
  std::array<std::array<float, 2>, 2> result = {0};
 
  // Start while loop for look for the zero
  auto current_step = first_guess;
  auto sign_alternate = 1;
  while (true)
  {
    auto current_delta = first_target->Eval(current_step) - second_target->Eval(current_step);
    current_step += starting_coarse;
  }
  return result;
}