UUBDownsampleFilter.h

Go to the documentation of this file.

#ifndef _sdet_UUBDownsampleFilter_h_
#define _sdet_UUBDownsampleFilter_h_

#include <utl/Trace.h>
#include <utl/TimeDistribution.h>
#include <utl/AugerUnits.h>
#include <utl/Math.h>


namespace sdet {

  namespace {

    // FIR coefficients from Dave Nitz's implementation in UUB firmware
    constexpr int kFirCoefficients[] = { 5, 0, 12, 22, 0, -61, -96, 0, 256, 551, 681, 551, 256, 0, -96, -61, 0, 22, 12, 0, 5 };
    // FIR normalization
    // in firmware an 11-bit right shift is used instead, ie 2048
    constexpr int kFirNormalizationBitShift = 11;
    //const int kFirNormalization = 2059;  // true norm of FIR
    //constexpr double kFirNormalization = (1 << kFirNormalizationBitShift);  // actually used
    constexpr double kUbSampling = 25*utl::nanosecond;
    constexpr int kADCSaturation = 4095;


    // enforce ADC saturation, both ways
    inline
    constexpr
    int
    Clip(const int i)
    {
      return std::max(0, std::min(i, kADCSaturation));
    }
 
  }


  // phase can be one only of { 0, 1, 2 }
  inline
  utl::TraceI
  UUBDownsampleFilter(const utl::TraceI& trace, const int phase = 1)
  {
    // input trace is assumed to have 8.333ns binning
    const int n = trace.GetSize();
    if (!n)
      return utl::TraceI(0, kUbSampling);
    const int m = utl::Length(kFirCoefficients);
    const int m2 = m / 2;
    std::vector<int> t;
    t.reserve(n + 2*m2);
    // pad front with the first trace values, but backwards
    for (int i = m2; i; --i)
      t.push_back(trace[i]);
    // copy the whole trace
    for (int i = 0; i < n; ++i)
      t.push_back(trace[i]);
    // pad back with the last trace values, but backwards
    for (int i = 1; i <= m2; ++i)
      t.push_back(trace[n-1-i]);
    const int n3 = n / 3;
    utl::TraceI res(n3, kUbSampling);
    for (int k = 0; k < n3; ++k) {
      auto& v = res[k];
      const int i = 3*k + phase;
      for (int j = 0; j < m; ++j)
        v += t[i + j] * kFirCoefficients[j];
      v >>= kFirNormalizationBitShift;
      v = Clip(v);
    }
    return res;
  }


  // phase can be one only of { 0, 1, 2 }
  inline
  utl::TimeDistributionI
  UUBDownsampleFilter(const utl::TimeDistributionI& trace, const int phase = 1)
  {
    const int n = trace.GetNumSlots();
    if (!n)
      return utl::TimeDistributionI(kUbSampling);
    const int m = utl::Length(kFirCoefficients);
    const int m2 = m / 2;
    std::vector<int> t;
    t.reserve(n + 2*m2);
    // pad front with first values, but backwards
    const int start = trace.GetStart();
    for (int i = m2; i; --i)
      t.push_back(trace.At(start + i));
    // copy the whole trace
    const int stop = trace.GetStop();
    for (int i = start; i <= stop; ++i)
      t.push_back(trace.At(i));
    // pad back with last values, but backwards
    for (int i = 1; i <= m2; ++i)
      t.push_back(trace.At(stop - i));
    utl::TimeDistributionI res(kUbSampling);
    const int startIndex = utl::FloorDiv(start+2, 3);
    const int stopIndex = utl::FloorDiv(stop-2, 3);
    for (int k = startIndex; k < stopIndex; ++k) {
      auto& v = res[k];
      const int i = 3*k + phase - start;
      for (int j = 0; j < m; ++j)
        v += t[i + j] * kFirCoefficients[j];
      v >>= kFirNormalizationBitShift;
      v = Clip(v);
    }
    return res;
  }


  // helpers for vectors of traces

  template<class Trace>
  inline
  std::vector<Trace>
  UUBDownsampleFilter(const std::vector<Trace>& traces)
  {
    std::vector<Trace> r;
    r.reserve(traces.size());
    for (const auto& t : traces)
      r.emplace_back(UUBDownsampleFilter(t));
    return r;
  }


  template<class Trace>
  inline
  std::vector<Trace>
  UUBDownsampleFilter(const std::vector<const Trace*>& traces)
  {
    std::vector<Trace> r;
    r.reserve(traces.size());
    for (const auto t : traces)
      r.emplace_back(UUBDownsampleFilter(*t));
    return r;
  }

}


#endif