TraceAlgorithm.cc

Go to the documentation of this file.

#include <utl/Trace.h>
#include <utl/TraceAlgorithm.h>
#include <utl/ErrorLogger.h>
#include <utl/AugerUnits.h>
#include <utl/Deprecator.h>
#include <cmath>
#include <vector>
#include <algorithm>
#include <numeric>
#include <sstream>

using namespace utl;
using namespace std;


template<typename T>
double
TraceAlgorithm::Min(const Trace<T>& trace,
                    const unsigned int bin1, const unsigned int bin2)
{
  if (!TraceAlgorithm::BoundsOk(trace, bin1, bin2))
    throw InvalidTraceBoundException();

  return *min_element(&trace[bin1], &trace[bin2] + 1);
}


template<typename T>
double
TraceAlgorithm::Max(const Trace<T>& trace,
                    const unsigned int bin1, const unsigned int bin2)
{
  if (!TraceAlgorithm::BoundsOk(trace, bin1, bin2))
    throw InvalidTraceBoundException();

  return *max_element(&trace[bin1], &trace[bin2] + 1);
}


template<typename T>
double
TraceAlgorithm::Mean(const Trace<T>& trace,
                     const unsigned int bin1, const unsigned int bin2)
{
  if (!TraceAlgorithm::BoundsOk(trace, bin1, bin2))
    throw InvalidTraceBoundException();

  return Sum(trace, bin1, bin2) / (bin2 + 1 - bin1);
}


template<typename T>
inline
double
TraceAlgorithm::Mean(const vector<T>& v)
{
  return accumulate(v.begin(), v.end(), T(0)) / double(v.size());
}


template<typename T>
double
TraceAlgorithm::Median(const Trace<T>& trace,
                       const unsigned int bin1, const unsigned int bin2,
                       const unsigned int sortAlgorithmLimit)
{
  if (!TraceAlgorithm::BoundsOk(trace, bin1, bin2))
    throw InvalidTraceBoundException();

  // work only on trace between bin1 and bin2
  vector<T> workingSet(&trace[bin1], &trace[bin2+1]);

  // Calculate median by sorting if the number of elements is small enough
  if (bin2 - bin1 <= sortAlgorithmLimit) {
    const size_t half = workingSet.size() / 2;
    // partially sort trace up to one element furhter than the middle
    partial_sort(workingSet.begin(), workingSet.begin() + half + 1, workingSet.end());

    // for odd sizes return the element in the middle
    // for even sizes return the mean of the two elements in the middle
    if (workingSet.size() & 1)
      return workingSet[half];
    else
      return 0.5*(workingSet[half - 1] + workingSet[half]);
  }

  // calculate starting values
  double mean = Mean(workingSet);
  // delta is the number of elements lower than the mean
  // which have to be eliminated
  double delta = 0.5 * workingSet.size();

  double median = 0;

  vector<T> lowerSet;
  vector<T> upperSet;
  do {
    // devide trace into subset of values lower and higher than the mean
    for (const auto& value : workingSet) {
      if (value < mean)
        lowerSet.push_back(value);
      else
        upperSet.push_back(value);
    }

    // choose the bigger set for the next round and look if the median is already found
    if (delta < lowerSet.size()) {
      // the elements of the upper set must be larger than the median
      // and are thus not considered anymore
      swap(workingSet, lowerSet);
      mean = Mean(workingSet);
      // Median is found if mean value is the maximum
      if (mean == *max_element(workingSet.begin(), workingSet.end())) {
        median = mean; // as mean == max
        break;
      }
    } else {
      // the elements of the lower set must be smaller than the median
      // and are thus eliminated for further searching
      delta -= lowerSet.size();
      swap(workingSet, upperSet);
      mean = Mean(workingSet);

      // Median is found if no values are left to eliminate
      // in this case, delta is 1/2 or 0
      // if delta = 0, the trace length was even and the median is taken
      // as the mean of the two median values.
      if (delta < 0.1) {          // if delta == 0
        median = 0.5*(*min_element(workingSet.begin(), workingSet.end()) +
                      *max_element(lowerSet.begin(), lowerSet.end()));
        break;
      }
      // Median is also found if mean value is the minimum or delta is 1/2
      if (mean == *min_element(workingSet.begin(), workingSet.end()) ||
          delta < 0.9) {     // if delta == 1/2
        median = *min_element(workingSet.begin(), workingSet.end());
        break;
      }
    }

    lowerSet.clear();
    upperSet.clear();

  } while (delta > 0);  // should not exit here, but when median is found

  return median;
}


template<typename T>
double
TraceAlgorithm::RMS(const Trace<T>& trace,
                    const unsigned int bin1, const unsigned int bin2)
{
  utl::Deprecator::GetInstance().Deprecated(
    "TraceAlgorithm::RMS()",
    "v2r8",
    "use TraceAlgorithm::RootMeanSquare() or TraceAlgorithm::StandardDeviation() "
    "depending on what you actually want"
  );

  if (!TraceAlgorithm::BoundsOk(trace, bin1, bin2))
    throw InvalidTraceBoundException();

  const double mean = Mean(trace, bin1, bin2);
  double squares = 0;
  for (unsigned int i = bin1; i <= bin2; ++i) {
    const double diff = trace[i] - mean;
    squares += diff * diff;
  }

  return sqrt(squares / (bin2 - bin1 + 1U));
}


template<typename T>
double
TraceAlgorithm::RootMeanSquare(const Trace<T>& trace,
                               const unsigned int bin1, const unsigned int bin2)
{
  if (!TraceAlgorithm::BoundsOk(trace, bin1, bin2))
    throw InvalidTraceBoundException();

  double squares = 0;
  for (unsigned int i = bin1; i <= bin2; ++i) {
    squares += trace[i]*trace[i];
  }

  return sqrt(squares / (bin2 - bin1 + 1U));
}


template<typename T>
double
TraceAlgorithm::StandardDeviation(const Trace<T>& trace,
                    const unsigned int bin1, const unsigned int bin2)
{
  if (!TraceAlgorithm::BoundsOk(trace, bin1, bin2))
    throw InvalidTraceBoundException();

  const double mean = Mean(trace, bin1, bin2);
  double squares = 0;
  for (unsigned int i = bin1; i <= bin2; ++i) {
    const double diff = trace[i] - mean;
    squares += diff * diff;
  }

  return sqrt(squares / (bin2 - bin1 + 1U));
}


template<typename T>
double
TraceAlgorithm::Sum(const Trace<T>& trace,
                    const unsigned int bin1, const unsigned int bin2)
{
  if (!TraceAlgorithm::BoundsOk(trace, bin1, bin2))
    throw InvalidTraceBoundException();

  return accumulate(&trace[bin1], &trace[bin2] + 1, T(0));
}


template<typename T>
double
TraceAlgorithm::Centroid(const Trace<T>& trace,
                         const unsigned int bin1, const unsigned int bin2)
{
  if (!TraceAlgorithm::BoundsOk(trace, bin1, bin2))
    throw InvalidTraceBoundException();

  double sum = 0;
  double moment = 0;

  for (unsigned int i = bin1; i <= bin2; ++i) {
    const double value = trace[i];
    sum += value;
    moment += value * (i + 0.5);
  }

  return moment / sum;
}


template<typename T>
bool
TraceAlgorithm::BoundsOk(const Trace<T>& trace,
                         const unsigned int bin1, const unsigned int bin2)
{
  if (bin2 >= trace.GetSize()) {
    ostringstream errorMsg;
    errorMsg << " Trace stops at " << trace.GetStop()
             << " while evaluation requested between " << bin1
             << " and " << bin2;
    ERROR(errorMsg);
    return false;
  }

  if (bin1 > bin2) {
    ostringstream errorMsg;
    errorMsg << "second bin specified can not be smaller than first: bin1 = "
             << bin1 << ", bin2 = " << bin2;
    ERROR(errorMsg);
    return false;
  }

  return true;
}


template<typename T>
double
TraceAlgorithm::ShapeParameter(const Trace<T>& trace,
                               const unsigned int bin1, const unsigned int bin2)
{
  // Definition see e.g. GAP Note 2003-076 by J. Cronin

  if (!TraceAlgorithm::BoundsOk(trace, bin1, bin2))
    throw InvalidTraceBoundException();

  const double binSize = trace.GetBinning();
  const double shapeParameterTime = 600 * nanosecond;
  const unsigned int shapeParameterBin =
    bin1 + (unsigned int)(shapeParameterTime / binSize);

  if (shapeParameterBin < bin2) {
    const double sumLate = Sum(trace, shapeParameterBin, bin2);
    const double sumEarly = Sum(trace, bin1, shapeParameterBin);
    return sumLate / sumEarly;
  }

  return 0;
}


template<typename T>
double
TraceAlgorithm::TimeAtRelativeSignalX(const Trace<T>& trace,
                                      const unsigned int bin1, const unsigned int bin2,
                                      const double percent)
{
  if (!TraceAlgorithm::BoundsOk(trace, bin1, bin2))
    throw InvalidTraceBoundException();

  const double sum = Sum(trace, bin1, bin2);
  if (sum < 0) {
    cerr << "Pmt trace seems wrong, signal negative (" << sum
         << ")..." << endl;
    return 0;
  }

  const double signal = sum * percent / 100.;

  double sumContent = 0;
  double sumContentPrev = 0;
  unsigned int i;
  for (i = bin1; i <= bin2; ++i) {
    sumContentPrev = sumContent;
    sumContent += trace[i];
    if (sumContent >= signal)
      break;
  }

  return ((signal - sumContent) / (sumContent - sumContentPrev) + i + 0.5) *
    trace.GetBinning();
}


template<typename T>
double
TraceAlgorithm::RiseTime(const Trace<T>& trace,
                         const unsigned int bin1, const unsigned int bin2)
{
  if (!TraceAlgorithm::BoundsOk(trace, bin1, bin2))
    throw InvalidTraceBoundException();

  return
    TimeAtRelativeSignalX(trace, bin1, bin2, 50.) -
    TimeAtRelativeSignalX(trace, bin1, bin2, 10.);
}


template<typename T>
double
TraceAlgorithm::FallTime(const Trace<T>& trace,
                         const unsigned int bin1, const unsigned int bin2)
{
  if (!TraceAlgorithm::BoundsOk(trace, bin1, bin2))
    throw InvalidTraceBoundException();

  return
    TimeAtRelativeSignalX(trace, bin1, bin2, 90.) -
    TimeAtRelativeSignalX(trace, bin1, bin2, 50.);
}


// explicit instantiations

template
double
TraceAlgorithm::Min(const Trace<int>& trace,
                    const unsigned int bin1, const unsigned int bin2);

template
double
TraceAlgorithm::Max(const Trace<int>& trace,
                    const unsigned int bin1, const unsigned int bin2);

template
double
TraceAlgorithm::Mean(const Trace<int>& trace,
                     const unsigned int bin1, const unsigned int bin2);

template
double
TraceAlgorithm::Median(const Trace<int>& trace,
                       const unsigned int bin1, const unsigned int bin2,
                       const unsigned int SortAlgorithmLimit);

template
double
TraceAlgorithm::RMS(const Trace<int>& trace,
                    const unsigned int bin1, const unsigned int bin2);

template
double
TraceAlgorithm::RootMeanSquare(const Trace<int>& trace,
                               const unsigned int bin1, const unsigned int bin2);

template
double
TraceAlgorithm::StandardDeviation(const Trace<int>& trace,
                                  const unsigned int bin1, const unsigned int bin2);

template
double
TraceAlgorithm::Centroid(const Trace<int>& trace,
                         const unsigned int bin1, const unsigned int bin2);

template
double
TraceAlgorithm::Sum(const Trace<int>& trace,
                    const unsigned int bin1, const unsigned int bin2);

template
double
TraceAlgorithm::ShapeParameter(const Trace<int>& trace,
                               const unsigned int bin1, const unsigned int bin2);

template
double
TraceAlgorithm::TimeAtRelativeSignalX(const Trace<int>& trace,
                                      const unsigned int bin1, const unsigned int bin2,
                                      const double percent);

template
double
TraceAlgorithm::RiseTime(const Trace<int>& trace,
                         const unsigned int bin1, const unsigned int bin2);

template
double
TraceAlgorithm::FallTime(const Trace<int>& trace,
                         const unsigned int bin1, const unsigned int bin2);

template
double
TraceAlgorithm::Min(const Trace<double>& trace,
                    const unsigned int bin1, const unsigned int bin2);

template
double
TraceAlgorithm::Max(const Trace<double>& trace,
                    const unsigned int bin1, const unsigned int bin2);

template
double
TraceAlgorithm::Mean(const Trace<double>& trace,
                     const unsigned int bin1, const unsigned int bin2);

template
double
TraceAlgorithm::Median(const Trace<double>& trace,
                       const unsigned int bin1, const unsigned int bin2,
                       const unsigned int sortAlgorithmLimit);

template
double
TraceAlgorithm::RMS(const Trace<double>& trace,
                    const unsigned int bin1, const unsigned int bin2);

template
double
TraceAlgorithm::RootMeanSquare(const Trace<double>& trace,
                               const unsigned int bin1, const unsigned int bin2);

template
double
TraceAlgorithm::StandardDeviation(const Trace<double>& trace,
                                  const unsigned int bin1, const unsigned int bin2);

template
double
TraceAlgorithm::Centroid(const Trace<double>& trace,
                         const unsigned int bin1, const unsigned int bin2);

template
double
TraceAlgorithm::Sum(const Trace<double>& trace,
                    const unsigned int bin1, const unsigned int bin2);

template
double
TraceAlgorithm::ShapeParameter(const Trace<double>& trace,
                               const unsigned int bin1, const unsigned int bin2);

template
double
TraceAlgorithm::TimeAtRelativeSignalX(const Trace<double>& trace,
                                      const unsigned int bin1, const unsigned int bin2,
                                      const double percent);

template
double
TraceAlgorithm::RiseTime(const Trace<double>& trace,
                         const unsigned int bin1, const unsigned int bin2);

template
double
TraceAlgorithm::FallTime(const Trace<double>& trace,
                         const unsigned int bin1, const unsigned int bin2);