testStationTriggerAlgorithm.cc

Go to the documentation of this file.

#include <tst/Verify.h>
#include <utl/AugerUnits.h>
#include <utl/TimeDistribution.h>
#include <utl/Trace.h>
#include <sdet/StationTriggerAlgorithm.h>
#include <cppunit/extensions/HelperMacros.h>
#include <vector>
#include <bitset>

using namespace std;
using namespace sdet;
using namespace utl;
using namespace tst;
using namespace sevt;

#define ASSERT_EQUAL(x, y) CPPUNIT_ASSERT(Verify<Equal>(x, y))


ostream&
operator<<(ostream& os, const StationTriggerInfo& info)
{
  return
    os << "("
          "PLD=" << bitset<16>(info.GetPLDBits()) << ", "
          "T2Th=" << info.IsT2Threshold() << ", "
          "start=" << info.GetStart() << ", "
          "latch=" << info.GetLatch() << ", "
          "stop=" << info.GetStop() << ", "
          "length=" << info.GetStop() - info.GetStart()
       << ')';
}


class TestStationTriggerAlgorithm : public CppUnit::TestFixture {

  CPPUNIT_TEST_SUITE(TestStationTriggerAlgorithm);
  CPPUNIT_TEST(TestOneT1ThresholdOnTimeDistributionI);
  CPPUNIT_TEST(TestOneT2ThresholdOnTimeDistributionI);
  CPPUNIT_TEST(TestManyT1ThresholdOnTimeDistributionI);
  //CPPUNIT_TEST(TestOneTOTOnTimeDistributionI);
  //CPPUNIT_TEST(TestManyTOTOnTimeDistributionI);
  //CPPUNIT_TEST(TestT1ThresholdAndTOTOnTimeDistributionI);
  //CPPUNIT_TEST(TestOneTOTdOnTimeDistributionI);
  //CPPUNIT_TEST(TestOneTOTdOnTimeDistributionIWithBaseline);
  //CPPUNIT_TEST(TestT1ThresholdOnTraceDWithBaseline);
  CPPUNIT_TEST_SUITE_END();

private:
  const double fUBSampling = 25*ns;
  const double fT1Threshold = 1.7;
  const double fTOTThreshold = 0.2;
  const double fTOTdThreshold = 0.2;
  const double fT2Threshold = 3.4;
  const double fTOTdDecayInBins = 67 / fUBSampling;
  const int fLatchBin = 246;
  const int fTraceLength = 3*1024/4;
  static bool fgShowLegend;

public:
  TestStationTriggerAlgorithm()
  {
    if (!fgShowLegend)
      return;
    fgShowLegend = false;
    cerr <<
      "       +----------------- ePLDRandom"           "\n"
      "       |"                                       "\n"
      "       |   +------------- ePLDTOTC (MoPS)"      "\n"
      "       |   |+------------ ePLDTOTB (TOTd)"      "\n"
      "       |   ||+----------- ePLDTOTA (TOT)"       "\n"
      "       |   |||+---------- ePLDThreshold"        "\n"
      "       |   ||||"                                "\n"
      "       |   ||||   +------ ePLDLatchRandom"      "\n"
      "       |   ||||   |+----- ePLDLatchTOTC (MoPS)" "\n"
      "       |   ||||   ||+---- ePLDLatchTOTB (TOTd)" "\n"
      "       |   ||||   |||+--- ePLDLatchTOTA (TOT)"  "\n"
      "       |   ||||   ||||+-- ePLDLatchThreshold"   "\n"
      "       |   ||||   |||||"                        "\n"
      "  (PLD=f...ba98...43210,"                       "\n"
      ;
  }

  void setUp() override { }

  void tearDown() override { }

  void
  TestOneT1ThresholdOnTimeDistributionI()
  {
    const double vem = 60;
    TimeDistributionI trace1(fUBSampling);
    TimeDistributionI trace2(fUBSampling);
    TimeDistributionI trace3(fUBSampling);
    const int bin = -123;
    trace1[bin] = trace2[bin] = trace3[bin] = 2*vem;

    const int start = bin - fLatchBin + 1;
    vector<StationTriggerInfo> info;
    info.emplace_back(StationTriggerData::ePLDLatchThreshold, false, start, bin, start + fTraceLength);

    Check(vem, 0, -1000, 2000, trace1, trace2, trace3, info, "one t1th");
  }

  void
  TestOneT2ThresholdOnTimeDistributionI()
  {
    const int vem = 60;
    TimeDistributionI trace1(fUBSampling);
    TimeDistributionI trace2(fUBSampling);
    TimeDistributionI trace3(fUBSampling);
    const int bin = -123;
    trace1[bin] = trace2[bin] = trace3[bin] = 3.5*vem;

    const int start = bin - fLatchBin + 1;
    vector<StationTriggerInfo> info;
    info.emplace_back(StationTriggerData::ePLDLatchThreshold, true, start, bin, start + fTraceLength);

    Check(vem, 0, -1000, 2000, trace1, trace2, trace3, info, "one t2th");
  }

  void
  TestManyT1ThresholdOnTimeDistributionI()
  {
    const int vem = 60;
    TimeDistributionI trace1(fUBSampling);
    TimeDistributionI trace2(fUBSampling);
    TimeDistributionI trace3(fUBSampling);

    const int bin[] = { -1000+300, 0+300, 1000+300 };
    constexpr int nBits = Length(bin);
    for (int i = 0; i < nBits; ++i) {
      trace1[bin[i]] = trace2[bin[i]] = trace3[bin[i]] = 2*vem;
      trace1[bin[i]+1] = trace2[bin[i]+1] = trace3[bin[i]+1] = 1*vem;
    }

    vector<StationTriggerInfo> info;
    for (int i = 0; i < nBits; ++i) {
      const int start = bin[i] - fLatchBin + 1;
      info.emplace_back(StationTriggerData::ePLDLatchThreshold, false, start, bin[i], start + fTraceLength);
    }

    Check(vem, 0, -1000, 2000, trace1, trace2, trace3, info, "many t1th");
  }

  void
  TestT1ThresholdAndTOTOnTimeDistributionI()
  {
    const int vem = 60;
    TimeDistributionI trace1(fUBSampling);
    TimeDistributionI trace2(fUBSampling);
    TimeDistributionI trace3(fUBSampling);

    // this should latch into T1 threshold
    const int thBin = -666;
    trace1[thBin] = trace2[thBin] = trace3[thBin] = 2*vem;

    // this should give post-latch TOT
    const int totBin = thBin + 200;
    const int nCoinc = 15;
    const int stride = 5; // increased stride so TOTd does not trigger.
    for (int i = 0; i < nCoinc; ++i) {
      const int bin = totBin + stride*i;
      // only 2-fold is required for TOT
      switch (abs(i%4)) {
      case 0: trace1[bin] = trace2[bin] = 0.3*vem; break;
      case 1: trace1[bin] = trace3[bin] = 0.3*vem; break;
      case 2: trace2[bin] = trace3[bin] = 0.3*vem; break;
      case 3: trace1[bin] = trace2[bin] = trace3[bin] = 0.3*vem ; break;
      }
    }

    const int start = thBin - fLatchBin + 1;
    vector<StationTriggerInfo> info;
    info.emplace_back(
      StationTriggerData::ePLDLatchThreshold | StationTriggerData::ePLDTOTA, false,
      start, thBin, start + fTraceLength
    );

    Check(vem, 0, -1000, 2000, trace1, trace2, trace3, info, "one t1th and one tot");
  }

  void
  TestOneTOTOnTimeDistributionI()
  {
    const int vem = 60;
    TimeDistributionI trace1(fUBSampling);
    TimeDistributionI trace2(fUBSampling);
    TimeDistributionI trace3(fUBSampling);

    // this should latch TOT
    const int totBin = 666;
    const int nCoinc = 20;
    const int stride = 5; // increased stride so TOTd does not trigger.
    for (int i = 0; i < nCoinc; ++i) {
      const int bin = totBin + stride*i;
      // only 2-fold is required for TOT
      switch (abs(i%4)) {
      case 0: trace1[bin] = trace2[bin] = 0.3*vem; break;
      case 1: trace1[bin] = trace3[bin] = 0.3*vem; break;
      case 2: trace2[bin] = trace3[bin] = 0.3*vem; break;
      case 3: trace1[bin] = trace2[bin] = trace3[bin] = 0.4*vem ; break;
      }
    }

    // trigger should fire on 13th
    const int latch = totBin + (13-1)*stride;
    const int start = latch - fLatchBin + 1;
    vector<StationTriggerInfo> info;
    info.emplace_back(
      StationTriggerData::ePLDLatchTOTA | StationTriggerData::ePLDTOTA, false,
      start, latch, start + fTraceLength
    );

    Check(vem, 0, -1000, 2000, trace1, trace2, trace3, info, "one tot");
  }

  void
  TestManyTOTOnTimeDistributionI()
  {
    const int vem = 60;
    TimeDistributionI trace1(fUBSampling);
    TimeDistributionI trace2(fUBSampling);
    TimeDistributionI trace3(fUBSampling);

    const int first = -1000;
    const int last = 2000;
    const int stride = 5; // increased stride so TOTd does not trigger.
    // this should latch TOT
    for (int i = first; i < last; i += stride) {
      // only 2-fold is required for TOT
      switch (abs(i % 4)) {
      case 0: trace1[i] = trace2[i] = 0.3*vem; break;
      case 1: trace1[i] = trace3[i] = 0.3*vem; break;
      case 2: trace2[i] = trace3[i] = 0.3*vem; break;
      case 3: trace1[i] = trace2[i] = trace3[i] = 0.4*vem ; break;
      }
    }

    const int length = last - first;
    const int traceLength = fTraceLength - (fLatchBin-1);
    const int nT1 = length/traceLength + bool((length % traceLength)/(12*3 + 1));

    vector<StationTriggerInfo> info;
    for (int i = 0; i < nT1; ++i) {
      const int latch = first + traceLength*i + (13-1)*stride;
      const int start = latch - fLatchBin + 1;
      info.emplace_back(
        StationTriggerData::ePLDLatchTOTA | StationTriggerData::ePLDTOTA, false,
        start, latch, start + fTraceLength
      );
    }

    Check(vem, 0, first, last, trace1, trace2, trace3, info, "many tot");
  }

  void
  TestOneTOTdOnTimeDistributionI()
  {
    const int vem = 60;
    TimeDistributionI trace1(fUBSampling);
    TimeDistributionI trace2(fUBSampling);
    TimeDistributionI trace3(fUBSampling);

    // this should latch TOTd
    const int totBin = 666;
    const int nCoinc = 20;
    const int stride = 2;
    for (int i = 0; i < nCoinc; ++i) {
      const int bin = totBin + stride*i;
      switch (abs(i%7)) {
      case 0: trace1[bin] = 0.22*vem; break;
      case 1: trace2[bin] = 0.22*vem; break;
      case 2: trace3[bin] = 0.22*vem; break;
      case 3: trace1[bin] = trace2[bin] = 0.22*vem; break;
      case 4: trace1[bin] = trace3[bin] = 0.22*vem; break;
      case 5: trace2[bin] = trace3[bin] = 0.22*vem; break;
      case 6: trace1[bin] = trace2[bin] = trace3[bin] = 0.25*vem ; break;
      }
    }

    // trigger should fire on 18th if multiplicity is 2
    const int latch = totBin + (18-1)*stride;
    // trigger should fire on 19th if multiplicity is 3
    // const int latch = totBin + (19 - 1)*stride;
    const int start = latch - fLatchBin + 1;
    vector<StationTriggerInfo> info;
    info.emplace_back(
      StationTriggerData::ePLDLatchTOTB | StationTriggerData::ePLDTOTB, false,
      start, latch, start + fTraceLength
    );
    Check(vem, 0, -1000, 2000, trace1, trace2, trace3, info, "two totd");

    // convolve with exponential response
    double previous_1 = 0;
    double previous_2 = 0;
    double previous_3 = 0;
    const double rate = exp(-1/fTOTdDecayInBins);
    for (int bin = totBin; bin < totBin + nCoinc*stride + 10; ++bin) {
      trace1[bin] += previous_1 * rate;
      trace2[bin] += previous_2 * rate;
      trace3[bin] += previous_3 * rate;
      previous_1 = trace1[bin];
      previous_2 = trace2[bin];
      previous_3 = trace3[bin];
    }

    info.clear();
    info.emplace_back(
      StationTriggerData::ePLDLatchTOTB | StationTriggerData::ePLDTOTB | StationTriggerData::ePLDTOTA, false,
      start, latch, start + fTraceLength
    );
    Check(vem, 0, -1000, 2000, trace1, trace2, trace3, info, "???");
  }

  void
  TestOneTOTdOnTimeDistributionIWithBaseline()
  {
    const int vem = 60;
    TimeDistributionI trace1(fUBSampling);
    TimeDistributionI trace2(fUBSampling);
    TimeDistributionI trace3(fUBSampling);

    // this should latch TOTd
    const int totBin = 666;
    const int nCoinc = 20;
    const int stride = 2;
    for (int i = 0; i < nCoinc; ++i) {
      const int bin = totBin + stride*i;
      switch (abs(i%7)) {
      case 0: trace1[bin] = 0.22*vem; break;
      case 1: trace2[bin] = 0.22*vem; break;
      case 2: trace3[bin] = 0.22*vem; break;
      case 3: trace1[bin] = trace2[bin] = 0.22*vem; break;
      case 4: trace1[bin] = trace3[bin] = 0.22*vem; break;
      case 5: trace2[bin] = trace3[bin] = 0.22*vem; break;
      case 6: trace1[bin] = trace2[bin] = trace3[bin] = 0.25*vem ; break;
      }
    }

    // convolve with exponential response
    double previous_1 = 0;
    double previous_2 = 0;
    double previous_3 = 0;
    const double rate = exp(-1/fTOTdDecayInBins);
    for (int bin = totBin; bin < totBin + nCoinc*stride + 10; ++bin) {
      trace1[bin] += previous_1 * rate;
      trace2[bin] += previous_2 * rate;
      trace3[bin] += previous_3 * rate;
      previous_1 = trace1[bin];
      previous_2 = trace2[bin];
      previous_3 = trace3[bin];
    }

    // add baseline
    const double baseline = 50;
    for (int bin = -1000; bin <= 2000; ++bin) {
      trace1[bin] += baseline;
      trace2[bin] += baseline;
      trace3[bin] += baseline;
    }

    // trigger should fire on 18th if multiplicity is 2
    const int latch = totBin + (18-1)*stride;
    // trigger should fire on 19th if multiplicity is 3
    // const int latch = totBin + (19 - 1)*stride;
    const int start = latch - fLatchBin + 1;

    vector<StationTriggerInfo> info;
    info.emplace_back(
      StationTriggerData::ePLDLatchTOTB | StationTriggerData::ePLDTOTB | StationTriggerData::ePLDTOTA, false,
      start, latch, start + fTraceLength
    );
    Check(vem, baseline, -1000, 2000, trace1, trace2, trace3, info, "one tot (with nonzero baseline)");
  }

  /*void
  TestT1ThresholdOnTraceDWithBaseline()
  {
    const double vem = 1;
    TraceD trace1(fTraceLength, fUBSampling);
    TraceD trace2(fTraceLength, fUBSampling);
    TraceD trace3(fTraceLength, fUBSampling);
    const int bin = fLatchBin - 1;
    trace1[bin] = trace2[bin] = trace3[bin] = 2*vem;

    const int start = 0;
    vector<StationTriggerInfo> info;
    info.emplace_back(
      StationTriggerData::ePLDLatchThreshold, false,
      start, bin, start + fTraceLength
    );

    // should produce one T1 threshold
    Check(vem, 0.1*vem, 0, fTraceLength, trace1, trace2, trace3, info);

    info.clear();

    // baseline is too high: no signal
    Check(vem, 1*vem, 0, fTraceLength, trace1, trace2, trace3, info, "t1th (on trace)");
  }*/

private:
  template<class Trace>
  void
  Check(const double vem, const double baseline,
        const int first, const int last,
        const Trace& t1, const Trace& t2, const Trace& t3,
        const vector<StationTriggerInfo>& shouldGet,
        const string& comment = "")
    const
  {
    if (!comment.empty())
      cerr << comment << '\n';

    using namespace sdet::Trigger;

    ThresholdUp t1th(int(round(baseline + fT1Threshold*vem)), 3, 3);  // multiplicity, n
    ThresholdUp t2th(int(round(baseline + fT2Threshold*vem)), 3, 3);  // multiplicity, n
    TimeOverThreshold tot(int(round(baseline + fTOTThreshold*vem)), 12, 120, 2, 3);  // occupancy, window, multiplicity, n
    DecayingIntegral integ(122, 0.5, 75., baseline, 3./4, 1./64, 3);  // window, delay, decay, n
    TimeOverThresholdDeconvolved totd(int(round(baseline + fTOTdThreshold*vem)), 45, 54,
                                      10, 120, 2, 3); // minCount, window, minIntegral, multiplicity, n
    MultiplicityOfPositiveSteps mops(4, 121, 3, 31, 3, // occupancy, window, minUp, maxUp, vetoOffset
                                     2, 3);  // minIntegral, multiplicity, n

    StationTriggerAlgorithm sTrig(t1th, t2th, tot, integ, totd, mops, fLatchBin, fTraceLength);

    const vector<const Trace*> traces = { &t1, &t2, &t3 };
    const vector<StationTriggerInfo> infos = sTrig.Run(first, last, traces);
    Match(infos, shouldGet);
    const auto stot = SimpleToT(fTOTThreshold*vem, baseline, 12, 120, 2, first, last, traces);
    cerr << "  SimpleToT " << stot.first << ' ' << stot.second << "\n\n";
  }

  static
  void
  Match(const vector<StationTriggerInfo>& infos1, const vector<StationTriggerInfo>& infos2)
  {
    if (infos1.size() != infos2.size()) {
      cerr << '\n';
      for (unsigned int i = 0; i < infos1.size(); ++i)
        cerr << "trigger info " << i << ": " << infos1[i] << '\n';
    }
    ASSERT_EQUAL(infos1.size(), infos2.size());
    for (unsigned int i = 0; i < infos1.size(); ++i) {
      const StationTriggerInfo& info1 = infos1[i];
      const StationTriggerInfo& info2 = infos2[i];
      cerr << "match " << i+1 << '/' << infos1.size() << ":\n  " << info1 << " =?=\n  " << info2;
      ASSERT_EQUAL(info1.GetPLDBits(), info2.GetPLDBits());
      ASSERT_EQUAL(info1.IsT1(), info2.IsT1());
      ASSERT_EQUAL(info1.IsT1Threshold(), info2.IsT1Threshold());
      ASSERT_EQUAL(info1.IsTimeOverThreshold(), info2.IsTimeOverThreshold());
      ASSERT_EQUAL(info1.IsT2(), info2.IsT2());
      ASSERT_EQUAL(info1.IsT2Threshold(), info2.IsT2Threshold());
      ASSERT_EQUAL(info1.GetTrigger(), info2.GetTrigger());
      ASSERT_EQUAL(info1.GetStart(), info2.GetStart());
      ASSERT_EQUAL(info1.GetLatch(), info2.GetLatch());
      ASSERT_EQUAL(info1.GetStop(), info2.GetStop());
      cerr << " OK\n";
    }
  }

  template<class Trace>
  static
  pair<bool, int>
  SimpleToT(const double totTh, const double baseline,
            const int occupancy, const int window, const int multiplicity,
            const int first, const int last,
            const vector<const Trace*>& traces)
  {
    const int maxi = max(last - window, first);
    for (int i = first; i <= maxi; ++i) {
      const int maxj = min(i + window, last);
      int trig = 0;
      for (int j = i; j <= maxj; ++j) {
        int m = 0;
        for (const auto t : traces) {
          const double s = (*t)[j] - baseline;
          if (s >= totTh)
            ++m;
        }
        if (m >= multiplicity) {
          ++trig;
          if (trig > occupancy)
            return make_pair(true, j);
        }
      }
    }
    return make_pair(false, 0);
  }

};

bool TestStationTriggerAlgorithm::fgShowLegend = true;


CPPUNIT_TEST_SUITE_REGISTRATION(TestStationTriggerAlgorithm);