StationTriggerAlgorithm.h

Go to the documentation of this file.

#ifndef _sdet_StationTriggerAlgorithm_h_
#define _sdet_StationTriggerAlgorithm_h_

#include <utl/Accumulator.h>
#include <sevt/StationTriggerData.h>
#include <sdet/StationTriggerInfo.h>
#include <sdet/UUBDownsampleFilter.h>
#include <utl/Math.h>
#include <vector>


namespace sdet {

  /*
    DV: In the firmware, various thresholds are keept as integers while
    baselines are implemented as fixed-point  numbers. The latter we here keep as
    doubles.  The input signal is always an integer representing the ADC counts.
  */

  namespace Trigger {

    // important note: "threshold" in all of the triggers is defined in such a
    // way that the threshold is considered passed when a trace value is
    // strictly larger (ie > and not >=) than the threshold.  the same holds
    // for other similar quantities, like minimal integral requirement

    class Threshold {
      // simple treshold trigger
    public:
      explicit
      Threshold(const int threshold,
                const unsigned int multiplicity, const unsigned int n) :
        fThreshold(n, threshold),
        fMultiplicity(multiplicity)
      { }

      template<class Array>
      explicit
      Threshold(const Array& thresholds,
                const unsigned int multiplicity, const unsigned int n) :
        fThreshold(&thresholds[0], &thresholds[n]),
        fMultiplicity(multiplicity)
      { }

      bool
      operator()(const int* const v)
        const
      {
        unsigned int coincidences = 0;
        for (unsigned int i = 0, n = fThreshold.size(); i < n; ++i)
          if (v[i] > fThreshold[i])
            ++coincidences;
        return coincidences >= fMultiplicity;
      }

    protected:
      std::vector<int> fThreshold;
      const unsigned int fMultiplicity;
    };


    class ThresholdUp : public Threshold {
      // the same as class Threshold but only on positive (up) slope
    public:
      explicit
      ThresholdUp(const int threshold,
                  const unsigned int multiplicity, const unsigned int n) :
        Threshold(threshold, multiplicity, n),
        fPrevious(n, false)
      { }

      template<class Array>
      explicit
      ThresholdUp(const Array& threshold,
                  const unsigned int multiplicity, const unsigned int n) :
        Threshold(threshold, multiplicity, n),
        fPrevious(n, false)
      { }

      bool
      operator()(const int* const v)
      {
        unsigned int coincidences = 0;
        for (unsigned int i = 0, n = Threshold::fThreshold.size(); i < n; ++i) {
          const auto t = (v[i] > Threshold::fThreshold[i]);
          auto& p = fPrevious[i];
          if (t && !p)
            ++coincidences;
          p = t;
        }
        return coincidences >= fMultiplicity;
      }

    private:
      std::vector<int> fPrevious;
    };


    template<typename T>
    class Fifo {
      // fifo to store adc history for various window operations
    public:
      Fifo() = default;

      explicit Fifo(const unsigned int length, const T& init = T()) : fFifo(length, init) { }

      void Assign(const unsigned int length, const T& init = T())
      { fFifo.assign(length, init); }

      T
      PushPop(const T& x)
      {
        fFifo.push_front(x);
        const auto y = fFifo.back();
        fFifo.pop_back();
        return y;
      }

      unsigned int GetLength() const { return fFifo.size(); }

    private:
      std::deque<T> fFifo;
    };


    class TimeOverThreshold {
      // time over threshold
    public:
      explicit
      TimeOverThreshold(const int threshold,
                        const int minOccupancy, const unsigned int window,
                        const unsigned int multiplicity, const unsigned int n) :
        fThreshold(threshold, multiplicity, n),
        fMinOccupancy(minOccupancy),
        fTriggersInWindow(window)
      { }

      template<class Vector>
      explicit
      TimeOverThreshold(const Vector& thresholds,
                        const int minOccupancy, const unsigned int window,
                        const unsigned int multiplicity, const unsigned int n) :
        fThreshold(thresholds, multiplicity, n),
        fMinOccupancy(minOccupancy),
        fTriggersInWindow(window)
      { }

      bool
      operator()(const int* const v)
      {
        // fifo
        const auto trigger = fThreshold(v);
        fCount += int(trigger) - int(fTriggersInWindow.PushPop(trigger));
        // note that >= would be more appropriate here but in the
        // firmware > is used with the fMinOccupancy definition
        return fCount > fMinOccupancy;
      }

      int GetCount() const { return fCount; }

    private:
      Threshold fThreshold;
      const int fMinOccupancy;
      int fCount = 0;
      Fifo<bool> fTriggersInWindow;
    };


    class DecayingIntegral {
      // leaky history integrator needed by ToTd and MoPS

      // note that the threshold for the integral is a double
    public:
      explicit
      DecayingIntegral(const unsigned int window,
                       const double decrement, const double threshold,
                       const double baseline,
                       const double baselineTolerance, const double baselineStep,
                       const unsigned int n) :
        fWindow(window),
        fIntegral(n),
        fIntegralQueue(n),
        fBaseline(n, baseline),
        fDecrement(decrement),
        fThreshold(threshold),
        fBaselineTolerance(baselineTolerance),
        fBaselineStep(baselineStep)
      { }

      template<class Array>
      explicit
      DecayingIntegral(const unsigned int window,
                       const double decrement, const double threshold,
                       const Array& baseline,
                       const double baselineTolerance, const double baselineStep,
                       const unsigned int n) :
        fWindow(window),
        fIntegral(n),
        fIntegralQueue(n),
        fBaseline(&baseline[0], &baseline[n]),
        fDecrement(decrement),
        fThreshold(threshold),
        fBaselineTolerance(baselineTolerance),
        fBaselineStep(baselineStep)
      { }

      void
      operator()(const int* const v)
      {
        const unsigned int n = fIntegral.size();
        // init on first call
        if (!fIntegralQueue.front().GetLength()) {
          for (unsigned int i = 0; i < n; ++i)
            fIntegralQueue[i].Assign(fWindow, 0);
        }
        for (unsigned int i = 0; i < n; ++i) {
          auto& baseline = fBaseline[i];
          const auto& value = v[i];
          UpdateBaseline(baseline, value);
          auto& currentIntegral = fIntegral[i];
          const auto windowIntegral = fIntegralQueue[i].PushPop(currentIntegral);
          currentIntegral += value - baseline - fDecrement;
          if (windowIntegral > fThreshold || currentIntegral < 0)
            currentIntegral = 0;
          const auto twiceThreshold = 2*fThreshold;
          if (currentIntegral > twiceThreshold)
            currentIntegral = twiceThreshold;
        }
      }

      bool HasTriggered(const unsigned int i) const { return fIntegral[i] > fThreshold; }

    private:
      void
      UpdateBaseline(double& baseline, const int v)
      {
        if (v > baseline) {
          if (v > baseline + fBaselineTolerance)
            baseline += 2*fBaselineStep;
          else
            baseline += fBaselineStep;
        } else if (v < baseline) {
          if (v < baseline - fBaselineTolerance)
            baseline -= 2*fBaselineStep;
          else
            baseline -= fBaselineStep;
        }
      }

      const unsigned int fWindow;
      std::vector<double> fIntegral;
      std::vector<Fifo<double>> fIntegralQueue;
      std::vector<double> fBaseline;
      const double fDecrement;
      const double fThreshold;
      const double fBaselineTolerance;
      const double fBaselineStep;
    };


    class TimeOverThresholdDeconvolved {
      // time over threshold with deconvolved trace

      // works on deconvolved trace; note that compared to the
      // TimeOverThreshold trigger, the coincidence and the multiplicity
      // conditions are swapped, ie n TimeOverThreshold triggers are run with
      // multiplicity=1 on n individual deconvolved traces and bound together
      // with the multiplicity condition given as the ctor parameter
    private:
      class Deconvolution {
      public:
        Deconvolution(const int decayFD, const int decayFN)
          : fFD(decayFD), fFN(decayFN) { }

        int
        operator()(const int current, const int previous)
          const
        {
          const int x = std::max(current*64 - previous*fFD, 0);
          constexpr int half = 1 << 9;
          return std::min((x*fFN + half) >> 10, 4096);
        }

      private:
        const int fFD;
        const int fFN;
      };

    public:
      explicit
      TimeOverThresholdDeconvolved(const int threshold,
                                   const int decayFD, const int decayFN,
                                   const unsigned int minCount, const unsigned int window,
                                   const unsigned int multiplicity, const unsigned int n) :
        fDeconvolve(decayFD, decayFN),
        fMultiplicity(multiplicity)
      {
        fTimeOverThreshold.reserve(n);
        for (unsigned int i = 0; i < n; ++i)
          fTimeOverThreshold.emplace_back(threshold, minCount, window, 1, 1);
      }

      template<class Array>
      explicit
      TimeOverThresholdDeconvolved(const Array& threshold,
                                   const int decayFD, const int decayFN,
                                   const unsigned int minCount, const unsigned int window,
                                   const unsigned int multiplicity, const unsigned int n) :
        fDeconvolve(decayFD, decayFN),
        fMultiplicity(multiplicity)
      {
        fTimeOverThreshold.reserve(n);
        for (unsigned int i = 0; i < n; ++i)
          fTimeOverThreshold.emplace_back(threshold[i], minCount, window, 1, 1);
      }

      bool
      operator()(const int* const v, const DecayingIntegral& integ)
      {
        const unsigned int n = fTimeOverThreshold.size();
        // init on first call
        if (fPrevious.empty())
          fPrevious.assign(v, v + n);
        unsigned int mult = 0;
        for (unsigned int i = 0; i < n; ++i) {
          const int value = v[i];
          const int decon = fDeconvolve(value, fPrevious[i]);
          fPrevious[i] = value;
          mult += (unsigned int)(fTimeOverThreshold[i](&decon) && integ.HasTriggered(i));
        }
        return mult >= fMultiplicity;
      }

    private:
      const Deconvolution fDeconvolve;
      std::vector<int> fPrevious;
      std::vector<TimeOverThreshold> fTimeOverThreshold;
      const unsigned int fMultiplicity;
    };


    class PositiveStepCounter {
      // mops helper class running on a single trace
    public:
      PositiveStepCounter(const int occupancy, const unsigned int window,
                          const int minUp, const int maxUp,
                          const int vetoOffset) :
        fOccupancy(occupancy),
        fWindow(window),
        fMinUp(minUp),
        fMaxUp(maxUp),
        fVetoOffset(vetoOffset)
      { }

      bool
      operator()(const int v)
      {
        if (fVeto) {
          fStart = v;
          --fVeto;
          return Update(v, false);
        } else if (v >= fPrevious) {
          // positive (or zero) increase
          return Update(v, false);
        } else {
          // negative drop
          const int jump = fPrevious - fStart;
          if (jump > 0)
            fVeto = std::max(0, IntLog2(jump) - fVetoOffset);
          fStart = v;
          return Update(v, fMinUp < jump && jump <= fMaxUp);
        }
      }

    private:
      bool
      Update(const int v, const bool isStep)
      {
        fPrevious = v;
        const auto prevCounts = fCounts;
        fCounts += int(isStep) - int(fWindow.PushPop(isStep));
        return prevCounts > fOccupancy;
      }

      static
      int
      IntLog2(unsigned int x)
      {
        if (x == 0)
          return -1;  // kind of -inf
#ifdef __GNUG__
        typedef std::numeric_limits<decltype(x)> nl;
        return nl::digits + nl::is_signed - 1 - __builtin_clz(x);
#else
        int r = 0;
        while (x >>= 1)
          ++r;
        return r;
#endif
      }

      const int fOccupancy;
      Fifo<bool> fWindow;
      const int fMinUp;
      const int fMaxUp;
      const int fVetoOffset;
      int fPrevious = 0;
      int fStart = 0;
      int fVeto = 1;
      int fCounts = 0;
    };


    class MultiplicityOfPositiveSteps {
      // multiplicity of positive steps
    public:
      MultiplicityOfPositiveSteps(const int occupancy, const unsigned int window,
                                  const int minUp, const int maxUp, const int vetoOffset,
                                  const unsigned int multiplicity,
                                  const unsigned int n) :  // =1 only for UB FADC traces
        fStepCounter(n, { occupancy, window, minUp, maxUp, vetoOffset }),
        fMultiplicity(multiplicity)
      { }

      bool
      operator()(const int* const v, const DecayingIntegral& integ)
      {
        unsigned int trig = 0;
        for (unsigned int i = 0, n = fStepCounter.size(); i < n; ++i)
          trig += (unsigned int)(fStepCounter[i](v[i]) && integ.HasTriggered(i));
        return trig >= fMultiplicity;
      }

    private:
      std::vector<PositiveStepCounter> fStepCounter;
      const unsigned int fMultiplicity;
    };

  }


  class StationTriggerAlgorithm {

  public:
    StationTriggerAlgorithm(const Trigger::ThresholdUp& t1th,
                            const Trigger::ThresholdUp& t2th,
                            const Trigger::TimeOverThreshold& tot,
                            const Trigger::DecayingIntegral& integ,
                            const Trigger::TimeOverThresholdDeconvolved& totd,
                            const Trigger::MultiplicityOfPositiveSteps& mops,
                            const int latchBin, const int traceLength) :
      fT1ThresholdTrigger(t1th),
      fT2ThresholdTrigger(t2th),
      fTOTTrigger(tot),
      fIntegral(integ),
      fTOTdTrigger(totd),
      fMOPSTrigger(mops),
      fLatchBin(latchBin),
      fTraceLength(traceLength)
    { }

    template<typename T, template<typename> class Trace>  // trace here is either Trace or TimeDistribution
    std::vector<StationTriggerInfo>
    Run(const int start, const int stop, const std::vector<const Trace<T>*>& traces)
    {
      const double binning = 25*utl::nanosecond;
      if (std::fabs(traces.front()->GetBinning() - binning)/binning < 1e-3)
        return RunUB(start, stop, fLatchBin, traces.size(), &traces.front(), fTraceLength);

      // UUB
      const double uubBinning = binning / 3;
      if (std::fabs(traces.front()->GetBinning() - uubBinning)/uubBinning > 1e-3) {
        std::ostringstream err;
        err << "The implemented trigger algorithms work correctly only for "
            << binning/utl::nanosecond << " ns or " << uubBinning/utl::nanosecond << " ns sampling rate "
               "while " << traces.front()->GetBinning()/utl::nanosecond << " ns was received";
        throw utl::DoesNotComputeException(err.str());
      }
      const auto downsampledTraces = UUBDownsampleFilter(traces);
      const auto newStart = utl::FloorDiv(start, 3);
      const auto newStop = utl::FloorDiv(stop, 3);
      auto triggerInfos = RunUB(newStart, newStop, fLatchBin/3, downsampledTraces, fTraceLength/3 + 1);
      // scale back to UUB
      for (auto& info : triggerInfos) {
        const int scaledStart = 3*info.GetStart();
        const int scaledLatch = 3*info.GetLatch();
        const int scaledStop = scaledStart + fTraceLength;
        info.SetStartLatchStop(scaledStart, scaledLatch, scaledStop);
      }
      return triggerInfos;
    }

  private:
    template<typename Trace>
    std::vector<StationTriggerInfo>
    RunUB(const int start, const int stop, const int latch,
          const std::vector<Trace>& traces, const int traceLength)
    {
      std::vector<const Trace*> tracePtrs;
      for (const auto& t : traces)
        tracePtrs.push_back(&t);
      return RunUB(start, stop, latch, tracePtrs.size(), &tracePtrs.front(), traceLength);
    }

    template<typename Trace>
    std::vector<StationTriggerInfo>
    RunUB(const int start, const int stop, const int latch,
          const int n, const Trace* const traces, const int traceLength)
    {
      std::vector<StationTriggerInfo> info;

      int i = start;
      do {
        int pld = sevt::StationTriggerData::ePLDNone;
        utl::Accumulator::Max<bool> isT2Th(false);
        // find latch triggers first
        for ( ; i < stop; ++i) {
          int values[n];
          for (int j = 0; j < n; ++j)
            values[j] = (*traces[j])[i];

          if (fT1ThresholdTrigger(values))
            pld |= sevt::StationTriggerData::ePLDLatchThreshold;
          isT2Th(fT2ThresholdTrigger(values));
          if (fTOTTrigger(values))
            pld |= sevt::StationTriggerData::ePLDLatchTOTA;
          fIntegral(values);
          if (fTOTdTrigger(values, fIntegral))
            pld |= sevt::StationTriggerData::ePLDLatchTOTB;
          if (fMOPSTrigger(values, fIntegral))
            pld |= sevt::StationTriggerData::ePLDLatchTOTC;
          // latch if any trigger found
          if (pld != sevt::StationTriggerData::ePLDNone)
            break;
        }

        // nothing found
        if (pld == sevt::StationTriggerData::ePLDNone)
          return info;

        // latch
        const int trigLatch = i;
        const int trigStart = trigLatch - latch + 1;
        const int trigStop = trigStart + traceLength;

        const int iMax = std::min(trigStop, stop);

        // any other triggers in the rest of the triggered trace?
        for (++i; i < iMax; ++i) {
          int values[n];
          for (int j = 0; j < n; ++j)
            values[j] = (*traces[j])[i];
          if (fT1ThresholdTrigger(values))
            pld |= sevt::StationTriggerData::ePLDThreshold;
          isT2Th(fT2ThresholdTrigger(values));
          if (fTOTTrigger(values))
            pld |= sevt::StationTriggerData::ePLDTOTA;
          fIntegral(values);
          if (fTOTdTrigger(values, fIntegral))
            pld |= sevt::StationTriggerData::ePLDTOTB;
          if (fMOPSTrigger(values, fIntegral))
            pld |= sevt::StationTriggerData::ePLDTOTC;
        }

        info.emplace_back(pld, isT2Th.GetMax(), trigStart, trigLatch, trigStop);

        // we do not simulate trigger dead-time although it definitely can be observed in the data

      } while (i < stop);

      return info;
    }

    Trigger::Threshold fT1ThresholdTrigger;
    Trigger::Threshold fT2ThresholdTrigger;
    Trigger::TimeOverThreshold fTOTTrigger;
    Trigger::DecayingIntegral fIntegral;
    Trigger::TimeOverThresholdDeconvolved fTOTdTrigger;
    Trigger::MultiplicityOfPositiveSteps fMOPSTrigger;
    int fLatchBin = 0;
    int fTraceLength = 0;

  };

}


#endif