utility.h

Go to the documentation of this file.

#pragma once
 
namespace utility
{
  //  Find first/lqast value above threshold
  //  --- negative_search triggers the search in negative values
  template <bool negative_search = false>
  std::pair<float, float>
  find_first_above(TGraph *gTarget, float threshold)
  {
    auto sign = 1 - negative_search * 2;
    for (int iPnt = 0; iPnt < gTarget->GetN(); iPnt++)
      if (sign * gTarget->GetPointY(iPnt) > sign * threshold)
        return {gTarget->GetPointX(iPnt), gTarget->GetPointY(iPnt)};
    return {-1., -1.};
  }
  template <bool negative_search = false>
  std::pair<float, float>
  find_last_above(TGraph *gTarget, float threshold)
  {
    auto sign = 1 - negative_search * 2;
    for (int iPnt = gTarget->GetN() - 1; iPnt >= 0; iPnt--)
      if (sign * gTarget->GetPointY(iPnt) > sign * threshold)
        return {gTarget->GetPointX(iPnt), gTarget->GetPointY(iPnt)};
    return {-1., -1.};
  }
 
  //  Fill a TProfile for the average of TGraphs
  void fill_profile(TProfile *hTarget, TGraph *gSource)
  {
    for (auto iPnt = 0; iPnt < gSource->GetN(); iPnt++)
    {
      hTarget->Fill(gSource->GetPointX(iPnt), gSource->GetPointY(iPnt));
    }
  }
  //  Fill a TH2F for the persistance of TGraphs
  void fill_persistance(TH2F *hTarget, TGraph *gSource)
  {
    for (auto iPnt = 0; iPnt < gSource->GetN(); iPnt++)
    {
      hTarget->Fill(gSource->GetPointX(iPnt), gSource->GetPointY(iPnt));
    }
  }
  //  Fill a TH2F for the persistance of TGraphs
  void fill_persistance(TH2F *hTarget, TF1 *fSource)
  {
    for (auto iBin = 1; iBin <= hTarget->GetNbinsX(); iBin++)
    {
      auto current_x = hTarget->GetXaxis()->GetBinCenter(iBin);
      hTarget->Fill(current_x, fSource->Eval(current_x));
    }
  }
  //  Fill TH1F from all points in a given interval
  //  --- Source is a TGraph
  void fill_from_interval(TH1F *hTarget, TGraph *gSource, float min, float max)
  {
    for (int iPnt = 0; iPnt < gSource->GetN(); iPnt++)
    {
      if (gSource->GetPointX(iPnt) < min)
        continue;
      if (gSource->GetPointX(iPnt) > max)
        break;
      hTarget->Fill(gSource->GetPointY(iPnt));
    }
  }
  //  --- Source is a TProfile
  void fill_from_interval(TH1F *hTarget, TProfile *gSource, float min, float max)
  {
    for (int iBin = 1; iBin <= gSource->GetNbinsX(); iBin++)
    {
      if (gSource->GetBinCenter(iBin) < min)
        continue;
      if (gSource->GetBinCenter(iBin) > max)
        break;
      hTarget->Fill(gSource->GetBinContent(iBin));
    }
  }
  //  Get Graph max
  template <bool negative_search = false>
  std::pair<float, float>
  get_graph_max(TGraph *gTarget, float min_point = -1., float max_point = -1.)
  {
    auto sign = 1 - negative_search * 2;
    auto result_x = gTarget->GetY()[0];
    auto result_y = gTarget->GetX()[0];
    for (int iPnt = 1; iPnt < gTarget->GetN(); iPnt++)
    {
      auto current_x = gTarget->GetX()[iPnt];
      if (((current_x < min_point) || (current_x > max_point)) && (min_point != max_point))
        continue;
      auto current_y = gTarget->GetY()[iPnt];
      if (sign * current_y > sign * result_y)
      {
        result_x = current_x;
        result_y = current_y;
      }
    }
    return {result_x, result_y};
  }
  //  Moving average of a graph
  template <int average_over_n_points = 2>
  TGraph *
  moving_average(TGraph *gTarget)
  {
    static_assert(average_over_n_points > 0, "averaging should be with at least one point");
    if (average_over_n_points == 1)
      return gTarget;
    auto gResult = new TGraph();
    gResult->SetName(gTarget->GetName());
    for (int iPnt = 0; iPnt < gTarget->GetN() - average_over_n_points; ++iPnt)
    {
      double new_x = 0.;
      double new_y = 0.;
      for (int jPnt = iPnt; jPnt < iPnt + average_over_n_points; ++jPnt)
      {
        new_x += gTarget->GetPointX(jPnt);
        new_y += gTarget->GetPointY(jPnt);
      }
      new_x /= (double)average_over_n_points;
      new_y /= (double)average_over_n_points;
      gResult->SetPoint(iPnt, new_x, new_y);
    }
    return gResult;
  }
  //  Max filter of a graph
  template <int average_over_n_points = 2>
  TGraph *
  maximum_filter(TGraph *gTarget)
  {
    static_assert(average_over_n_points > 0, "filtering should be with at least one point");
    if (average_over_n_points == 1)
      return gTarget;
    auto gResult = new TGraph();
    gResult->SetName(gTarget->GetName());
    auto new_points = 0;
    for (int iPnt = 0; iPnt < gTarget->GetN(); iPnt += average_over_n_points)
    {
      double new_x = gTarget->GetPointX(iPnt);
      double new_y = gTarget->GetPointY(iPnt);
      for (int jPnt = iPnt + 1; jPnt < min(iPnt + average_over_n_points,gTarget->GetN()); ++jPnt)
      {
        if (gTarget->GetPointY(jPnt) > new_y)
        {
          new_x = gTarget->GetPointX(jPnt);
          new_y = gTarget->GetPointY(jPnt);
        }
      }
      gResult->SetPoint(new_points, new_x, new_y);
      new_points++;
    }
    return gResult;
  }
  //  Gauss filter of a graph
  template <int window_of_n_points = 2>
  TGraph *
  gauss_filter(TGraph *gTarget, float alpha_parameter = 2.5)
  {
    //  Require the window of points to be at least 1
    static_assert(window_of_n_points > 0, "gaus filter should be with at least one point");
    //  >>> If it is one, return target itself
    if (window_of_n_points == 1)
      return gTarget;
    //  Build result graph
    auto gResult = new TGraph();
    gResult->SetName(Form("%s_gausfilter_n%i_a%i", gTarget->GetName(), window_of_n_points));
    //
    for (int iPnt = 0; iPnt < gTarget->GetN() - window_of_n_points; ++iPnt)
    {
      double new_x = 0.;
      double new_y = 0.;
      for (int jPnt = iPnt; jPnt < iPnt + window_of_n_points; ++jPnt)
      {
        new_x += gTarget->GetPointX(jPnt);
        new_y += gTarget->GetPointY(jPnt);
      }
      new_x /= (double)window_of_n_points;
      new_y /= (double)window_of_n_points;
      gResult->SetPoint(iPnt, new_x, new_y);
    }
    return gResult;
  }
  //  Derivative of graph
  TGraph *
  derivate(TGraph *gTarget, double sign = 1.)
  {
    auto gResult = new TGraph();
    gResult->SetName(gTarget->GetName());
    for (int iPnt = 0; iPnt < gTarget->GetN() - 1; ++iPnt)
    {
      auto current_x = gTarget->GetPointX(iPnt);
      auto current_y = gTarget->GetPointY(iPnt);
      auto next_x = gTarget->GetPointX(iPnt + 1);
      auto next_y = gTarget->GetPointY(iPnt + 1);
      gResult->SetPoint(iPnt, current_x * 0.5 + next_x * 0.5, (next_y - current_y) / (next_x - current_x));
    }
    return gResult;
  }
  //  From TProfile to TGraphErrors
  TGraph *
  TProfile_to_TGraphErorrs(TProfile *gSource)
  {
    auto result = new TGraphErrors();
    for (auto iBin = 1; iBin < gSource->GetNbinsX(); iBin++)
    {
      result->SetPoint(iBin - 1, gSource->GetBinCenter(iBin), gSource->GetBinContent(iBin));
      result->SetPointError(iBin - 1, gSource->GetBinWidth(iBin) * 0.5, gSource->GetBinError(iBin));
    }
    return result;
  }
  //  Measure simil-integral from waveform
  template <int power = 1,
            bool bin_width = false>
  float
  integrate(TGraph *gTarget, float min_point = -1., float max_point = -1.)
  {
    auto gResult = new TGraph();
    gResult->SetName(gTarget->GetName());
    auto bin_size = bin_width ? (gTarget->GetPointX(1) - gTarget->GetPointX(0)) : 1.;
    auto integral = 0.;
    for (int iPnt = 0; iPnt < gTarget->GetN() - 1; ++iPnt)
    {
      auto current_x = gTarget->GetX()[iPnt];
      if (((current_x < min_point) || (current_x > max_point)) && (min_point != max_point))
        continue;
      integral += bin_size * pow(gTarget->GetPointY(iPnt), power);
    }
    return integral;
  }
  //  Measure noisepower
  //  --- From TGraph
  template <int n_parameter = 2,
            int n_averaging = 10>
  std::pair<float, float>
  measure_noisepower(TGraph *gTarget)
  {
    //  Derivate the requeste number of time, smoothing over the given number of averging points
    static_assert(n_parameter > 0, "n parameter must be greater than 0");
    static_assert(n_averaging > 0, "averaging must be greater than 0");
    TGraph *used_derivative;
    used_derivative = derivate(moving_average<n_averaging>(gTarget));
    for (auto iTer = 1; iTer < n_parameter; iTer++)
      used_derivative = derivate(moving_average<n_averaging>(used_derivative));
    //
    return {integrate<2>(used_derivative, -2., 2.) - 0.04 * integrate<2>(used_derivative, -10., -110.), integrate(gTarget, -5., 95.) - integrate(gTarget, -10., -110.)};
  }
  //  --- From TProfile
  template <int n_parameter = 2,
            int n_averaging = 10>
  std::pair<float, float>
  measure_noisepower(TProfile *hTarget)
  {
    return measure_noisepower<n_parameter, n_averaging>(TProfile_to_TGraphErorrs(hTarget));
  }
 
  //  Peak finder
  //  Each element is a peak that has been found:
  //  [peak x-coordinate] [peak y-coordinate] [point n]
  /*
  template <int window_of_n_points = 100,
            int average_over_n_points = 100,
            bool negative_search = false>
  std::vector<std::array<float, 3>>
  find_peaks(TGraph *gTarget, float min_height = 10.)
  {
    //  Result
    std::vector<std::array<float, 3>> result;
    //  Check the window is at least 3
    static_assert(window_of_n_points > 2, "peak finder should have at least 3 points window");
    auto gNewTarget = moving_average<average_over_n_points>(gTarget);
    for (int iPnt = 0; iPnt < gNewTarget->GetN(); iPnt++)
    {
      auto current_x = gNewTarget->GetPointX(iPnt);
      auto current_y = gNewTarget->GetPointY(iPnt);
      auto previous_y = gNewTarget->GetPointY(max(0, iPnt - window_of_n_points));
      auto next_y = gNewTarget->GetPointY(min(gNewTarget->GetN(), iPnt + window_of_n_points));
      if ((current_y - previous_y) < min_height)
        continue;
      if ((current_y - next_y) < min_height)
        continue;
      result.push_back({static_cast<float>(current_x), static_cast<float>(current_y), static_cast<float>(iPnt)});
      iPnt += window_of_n_points;
    }
    TCanvas* c1 = new TCanvas();
    gNewTarget->Draw("ALPE");
    return result;
  }
  */
 
  template <int window_of_n_points = 10,
            int average_over_n_points = 20,
            bool negative_search = false>
  std::vector<std::array<float, 3>>
  find_peaks(TGraph *gTarget, float min_height = 10.)
  {
    //  Result
    std::vector<std::array<float, 3>> result;
    //  Check the window is at least 3
    static_assert(window_of_n_points > 2, "peak finder should have at least 3 points window");
    auto gNewTarget = maximum_filter<average_over_n_points>(gTarget);
    for (int iPnt = 0; iPnt < gNewTarget->GetN(); iPnt++)
    {
      auto current_x = gNewTarget->GetPointX(iPnt);
      auto current_y = gNewTarget->GetPointY(iPnt);
      auto previous_y = gNewTarget->GetPointY(max(0, iPnt - window_of_n_points));
      auto next_y = gNewTarget->GetPointY(min(gNewTarget->GetN(), iPnt + window_of_n_points));
      if ((current_y - previous_y) < min_height)
        continue;
      if ((current_y - next_y) < min_height)
        continue;
      result.push_back({static_cast<float>(current_x), static_cast<float>(current_y), static_cast<float>(iPnt)});
      iPnt += window_of_n_points;
    }
    TCanvas *c1 = new TCanvas();
    gNewTarget->Draw("ALPE");
    return result;
  }
 
} // namespace utility