SdEventSelector.cc

Go to the documentation of this file.

#include "SdEventSelector.h"

#include <fwk/CentralConfig.h>
#include <fwk/RunController.h>

#include <det/Detector.h>

#include <sdet/SDetector.h>
#include <sdet/Station.h>

#include <evt/Event.h>
#include <evt/ShowerRecData.h>
#include <evt/ShowerSRecData.h>

#include <sevt/SEvent.h>
#include <sevt/Header.h>
#include <sevt/EventTrigger.h>
#include <sevt/Station.h>
#include <sevt/StationTriggerData.h>
#include <sevt/StationSimData.h>
#include <sevt/StationRecData.h>
#include <sevt/PMT.h>
#include <sevt/PMTRecData.h>
#include <sevt/PMTCalibData.h>
#include <sevt/IsLightning.h>

#include <utl/Vector.h>
#include <utl/AxialVector.h>
#include <utl/GeometryUtilities.h>
#include <utl/PhysicalConstants.h>
#include <utl/ErrorLogger.h>
#include <utl/String.h>
#include <utl/Reader.h>
#include <utl/TabularStream.h>

#include <TMinuit.h>

#include <boost/range/iterator_range.hpp>

#include <set>

using namespace fwk;
using namespace evt;
using namespace sevt;
using namespace utl;
using namespace std;
using namespace SdEventSelectorOG;
using boost::make_iterator_range;


// globals

namespace SdEventSelectorOG {

  /*
    In triangular grid distances grow as following

                    1500 m
    shell fac       array
    0     0         0
    1     1         1500      C1
    2     1.73205   2598      hexagon skew diagonal
    3     2         3000      hexagon diameter
    4     2.64575   3969
    5     3         4500
    6     3.46410   5196
    7     3.60555   5408
  */

  //inlines

  // stl helpers
  template<typename Iterator>
  inline
  Iterator
  Next(Iterator it)
  {
    return ++it;
  }


  // distance helper
  inline
  double
  Distance(const sdet::Station& i, const sdet::Station& j)
  {
    return (i.GetPosition() - j.GetPosition()).GetMag();
  }


  inline
  bool
  IsInRange(const double x, const vector<double>& r)
  {
    return (r[0] < x && x < r[1]);
  }


  namespace SeedFit {

    double gScalarProductIJ;
    double gScalarProductJK;
    double gTimeI;
    double gTimeJ;

    void
    Chi2(int& /*nPar*/, double* const /*grad*/, double& value,
         double* const par, const int /*flag*/)
    {
      const double& phi = par[0];

      const double u = cos(phi);
      const double v = sin(phi);

      value =
        Sqr(u - gTimeI) +
        Sqr(u * gScalarProductIJ + v * gScalarProductJK - gTimeJ);
    }

  } // namespace SeedFit


  ShowerSRecData&
  GetSRecShower(evt::Event& event)
  {
    if (!event.HasRecShower())
      event.MakeRecShower();
    auto& shRec = event.GetRecShower();
    if (!shRec.HasSRecShower())
      shRec.MakeSRecShower();
    return shRec.GetSRecShower();
  }


  void
  MakeInfo(const vector<int>& ids, const string& message)
  {
    if (ids.empty())
      return;
    ostringstream info;
    info << "station" << String::Plural(ids) << ' ' << String::OfSortedIds(ids) << ' ' << message;
    INFO(info);
  }


  inline
  double
  SdEventSelector::CTimeDifferenceMargin(const Station& i, const Station& j)
    const
  {
    return kSpeedOfLight*(abs((i.GetRecData().GetSignalStartTime() -
                               j.GetRecData().GetSignalStartTime()).GetInterval()) -
                          fLightCompatibilityTolerance);
  }


  int
  SdEventSelector::StationSelection(sevt::SEvent& sEvent)
    const
  {
    vector<int> ssElectronicsType;
    vector<int> ssOffGrid;
    vector<int> ssDense;
    vector<int> ssLightning;
    vector<int> ssNoRecData;
    vector<int> ssTOTd;
    vector<int> ssMoPS;

    int nCandidates = 0;
    for (auto& s : sEvent.StationsRange()) {
      const int sId = s.GetId();
      const auto& dStation = det::Detector::GetInstance().GetSDetector().GetStation(s);

      if (fRejectUBStations && !s.IsUUB()) {
        ssElectronicsType.push_back(sId);
        s.SetRejected(StationConstants::eElectronicsType);
      }

      if (fRejectUUBStations && s.IsUUB()) {
        ssElectronicsType.push_back(sId);
        s.SetRejected(StationConstants::eElectronicsType);
      }

      if (!dStation.IsInGrid(sdet::SDetectorConstants::GridIndex(fGridType))) {
        ssOffGrid.push_back(sId);
        s.SetRejected(StationConstants::eOffGrid);
      }

      if (dStation.IsDense()) {
        ssDense.push_back(sId);
        s.SetRejected(StationConstants::eDenseArray);
      }

      if (s.IsCandidate()) {

        if (!s.IsUUB()) {

          /* For now we skip upgraded stations
             Stations with new electronics are wrongly flagged as lightning
             as the algorithm does not handle properly new traces (higher sampling rate and resolution in amplitude)
             TODO: update Lightning algorithm for UUB
          */

          if (IsLightning(s)) {
            ssLightning.push_back(sId);
            switch (fRejectLightning) {
            case 1:
              // flag station
              s.SetRejected(StationConstants::eLightning);
              break;
            case 2:
              // reject whole event
              WARNING("Lightning detected, whole event rejected.");
              return 0;
            }
          }
        }

        if (!s.HasRecData()) {
          ssNoRecData.push_back(sId);
          s.SetRejected(StationConstants::eNoRecData);
        }

        const auto& sTrig = s.GetTriggerData();
        const bool thtot = sTrig.IsThresholdOrTimeOverThreshold(); // old triggers

        if (!thtot) {
          const bool totd = sTrig.IsTimeOverThresholdDeconvoluted();
          const bool silenceTotd = fSilenceTOTdStations && totd;
          const bool mops = sTrig.IsMultiplicityOfPositiveSteps();
          const bool silenceMops = fSilenceMoPSStations && mops;
          if ((fSilenceTOTdStations && silenceMops) || (silenceTotd && !mops) || (silenceMops && !totd)) {
            if (!s.IsRejected())
              s.SetSilent();
            if (totd)
              ssTOTd.push_back(sId);
            if (mops)
              ssMoPS.push_back(sId);
          }
        }

        ++nCandidates;

      }
    }

    MakeInfo(ssOffGrid, "rejected as off-grid.");
    MakeInfo(ssDense, "rejected as dense.");
    MakeInfo(ssLightning, fRejectLightning ? "rejected as potential lightning." : "NOT rejected as lightning.");
    MakeInfo(ssElectronicsType, "rejected due to electronics type");
    MakeInfo(ssNoRecData, "without RecData rejected.");
    MakeInfo(ssTOTd, "set silent because only TOTd.");
    MakeInfo(ssMoPS, "set silent because only MoPS.");

    return nCandidates;
  }


  string
  GetT4TriggerName(const bool isT4,
                   const bool isFD,
                   const bool is3TOT,
                   const bool is3TOTd,
                   const bool is3MoPS,
                   const bool is3TOTMix,
                   const bool is4C1)
  {
    ostringstream trig;
    if (isT4) {
      if (isFD)
        trig << "FD";
      if (is3TOT)
        trig << (trig.str().length() ? "+" : "") << "3TOT";
      if (is3TOTd)
        trig << (trig.str().length() ? "+" : "") << "3TOTd";
      if (is3MoPS)
        trig << (trig.str().length() ? "+" : "") << "3MoPS";
      if (!is3TOT && !is3TOTd && !is3MoPS && is3TOTMix)
        trig << (trig.str().length() ? "+" : "") << "3TOTMix";
      if (is4C1)
        trig << (trig.str().length() ? "+" : "") << "4C1";
    } else
      trig << "None";
    return trig.str();
  }


  VModule::ResultFlag
  SdEventSelector::Init()
  {
    auto topB = CentralConfig::GetInstance()->GetTopBranch("SdEventSelector");

    topB.GetChild("MinNumberOfStationsAtBeginning").GetData(fMinNumberOfStationsAtBeginning);

    topB.GetChild("EnableBottomUpSelection").GetData(fEnableBottomUpSelection);

    topB.GetChild("EnableT4Trigger").GetData(fEnableT4Trigger);

    topB.GetChild("RejectNonT4Events").GetData(fRejectNonT4Events);

    topB.GetChild("EnableT5Trigger").GetData(fEnableT5Trigger);

    topB.GetChild("MinNumberOfActiveStations").GetData(fMinNumberOfActiveStations);

    topB.GetChild("RejectNonT5Events").GetData(fRejectNonT5Events);

    topB.GetChild("RejectLightning").GetData(fRejectLightning);
    switch (fRejectLightning) {
    case 0:
      INFO("NOT rejecting stations with potential lightning signal.");
      break;
    case 1:
      INFO("Rejecting individual stations with lightning signal.");
      break;
    case 2:
      INFO("Rejecting whole event when lightning detected.");
    }

    topB.GetChild("RejectUBStations").GetData(fRejectUBStations);
    topB.GetChild("RejectUUBStations").GetData(fRejectUUBStations);

    if (fRejectUBStations && fRejectUUBStations) {
      ERROR("Both UB and UUB stations configured for rejection. "
            "At least one type of electronics must remain selected.");
      return eFailure;
    }

    // bad periods
    const bool bad = topB.GetChild("EnableSdExcludedPeriods").Get<bool>();

    fBadPeriods.clear();

    if (bad) {
      topB.GetChild("RejectBadPeriods").GetData(fRejectBadPeriods);

      if (fRejectBadPeriods)
        INFO("Events in bad periods are rejected.");
      else
        INFO("Events in bad periods are tagged.");

      auto excludedB = CentralConfig::GetInstance()->GetTopBranch("SdExcludedPeriods");

      for (auto periodB = excludedB.GetFirstChild();
           periodB; periodB = periodB.GetNextSibling()) {
        const auto& att = periodB.GetAttributes();
        const int id = boost::lexical_cast<int>(att.find("id")->second);
        const TimeStamp begin = periodB.GetChild("begin").Get<TimeStamp>();
        const TimeStamp end = periodB.GetChild("end").Get<TimeStamp>();
        fBadPeriods.push_back(SdEventSelectorOG::BadPeriod(id, TimeRange(begin, end)));
      }

      if (fBadPeriods.begin() == fBadPeriods.end())
        WARNING("No excluded periods found.");
      else {
        sort(fBadPeriods.begin(), fBadPeriods.end());

        auto bPrevious = fBadPeriods.begin();
        const auto bEnd = fBadPeriods.end();
        if (bPrevious != bEnd) {
          auto bIt = bPrevious;
          for (++bIt; bIt != bEnd; ++bIt) {
            if (bPrevious->HasCommonTime(*bIt)) {
              ostringstream err;
              err << "Excluded period "
                     "id=" << bPrevious->GetId() << " "
                     "[" << bPrevious->GetStartTime() << " | " << bPrevious->GetStopTime() << "] "
                     "is in conflict with period "
                     "id=" << bIt->GetId() <<" "
                     "[" << bIt->GetStartTime() << " | " << bIt->GetStopTime() << "].";
              ERROR(err);
              return eFailure;
            }
            bPrevious = bIt;
          }
        }
      }
    }

    // selection parameters
    topB.GetChild("GridType").GetData(fGridType);
    switch (fGridType) {
    case sdet::SDetectorConstants::eStandard:
      INFO("Selecting stations from 1500 m array.");
      break;
    case sdet::SDetectorConstants::eInfill750:
      INFO("Selecting stations from 750 m array.");
      break;
    case sdet::SDetectorConstants::eInfill433:
      INFO("Selecting stations from 433 m array.");
      break;
    default:
      ERROR("Unknown GridType!");
      return eFailure;
    }

    topB.GetChild("LightCompatibilityTolerance").GetData(fLightCompatibilityTolerance);
    topB.GetChild("FirstCrownDistanceRange").GetData(fFirstCrownDistanceRange);
    topB.GetChild("SkewMaximumDistance").GetData(fSkewMaximumDistance);
    topB.GetChild("SeedPlaneFrontTimeResidualRangeToT").GetData(fSeedPlaneFrontTimeResidualRangeToT);
    topB.GetChild("SeedPlaneFrontTimeResidualRangeToTd").GetData(fSeedPlaneFrontTimeResidualRangeToTd);
    topB.GetChild("SeedPlaneFrontTimeResidualRangeMoPS").GetData(fSeedPlaneFrontTimeResidualRangeMoPS);
    topB.GetChild("SeedPlaneFrontTimeResidualRangeTh2").GetData(fSeedPlaneFrontTimeResidualRangeTh2);
    topB.GetChild("SeedPlaneFrontTimeResidualRangeTh1").GetData(fSeedPlaneFrontTimeResidualRangeTh1);
    topB.GetChild("LonelyIfNoneInDistance").GetData(fLonelyIfNoneInDistance);
    topB.GetChild("LonelyIfOneInDistance").GetData(fLonelyIfOneInDistance);

    topB.GetChild("SilenceTOTdStations").GetData(fSilenceTOTdStations);
    if (fSilenceTOTdStations)
      INFO("Silencing stations with ToTd-only trigger.");

    topB.GetChild("SilenceMoPSStations").GetData(fSilenceMoPSStations);
    if (fSilenceMoPSStations)
      INFO("Silencing stations with MoPS-only trigger.");

    const auto silenceNewB = topB.GetChild("SilenceTOTdMoPSTriggersBeforeDate");
    if (silenceNewB.GetChild("Enabled").Get<bool>()) {
      silenceNewB.GetChild("Date").GetData(fSilenceTOTdMoPSDate);
      ostringstream info;
      info << "Will silence stations with ToTd or MoPS only trigger before the date " << fSilenceTOTdMoPSDate;
      INFO(info);
    }

    return eSuccess;
  }


  VModule::ResultFlag
  SdEventSelector::Run(evt::Event& event)
  {
    auto& namedCounters = RunController::GetInstance().GetRunData().GetNamedCounters();

    if (!event.HasSEvent()) {
      ++namedCounters["SdEventSelector/NoSdEvent"];
      return eContinueLoop;
    }

    auto& sEvent = event.GetSEvent();

    // station selection
    const auto nStations = StationSelection(sEvent);
    if (nStations < fMinNumberOfStationsAtBeginning) {
      ostringstream info;
      info << "Not enough stations (" << nStations << " < " << fMinNumberOfStationsAtBeginning << ")!";
      INFO(info);
      ++namedCounters["SdEventSelector/NotEnoughStations"];
      return eContinueLoop;
    }

    fCandidateStations.clear();
    f3TOTCandidateStations.clear();
    f3TOTdCandidateStations.clear();
    f3MoPSCandidateStations.clear();
    f3TOTMixCandidateStations.clear();
    f4C1CandidateStations.clear();
    fSeedStations.clear();

    const bool isFD = sEvent.GetTrigger().IsFD();

    if (isFD)
      ++namedCounters["SdEventSelector/FdTrigger"];

    const auto& eventTime = sEvent.GetHeader().GetTime();
    if (fSilenceTOTdMoPSDate) {
      if (eventTime < fSilenceTOTdMoPSDate) {
        fSilenceTOTdStations = true;
        fSilenceMoPSStations = true;
      } else {
        fSilenceTOTdStations = false;
        fSilenceMoPSStations = false;
      }
    }

    auto& shSRec = GetSRecShower(event);
    shSRec.SetIsTOTdMoPSSilent(fSilenceTOTdStations && fSilenceMoPSStations);

    const auto& sDet = det::Detector::GetInstance().GetSDetector();

    // prepare station lists
    for (auto& s : sEvent.CandidateStationsRange()) {
      const auto& sTrig = s.GetTriggerData();
      const auto& ds = sDet.GetStation(s);
      const StationPair sp(&s, &ds);
      fCandidateStations.push_back(sp);
      if (sTrig.IsT2() || isFD) {
        f4C1CandidateStations.push_back(sp);
        if (sTrig.IsTimeOverThreshold())
          f3TOTCandidateStations.push_back(sp);
        if (sTrig.IsTimeOverThresholdDeconvoluted() && !fSilenceTOTdStations)
          f3TOTdCandidateStations.push_back(sp);
        if (sTrig.IsMultiplicityOfPositiveSteps() && !fSilenceMoPSStations)
          f3MoPSCandidateStations.push_back(sp);
        if (sTrig.IsTimeOverThreshold() ||
            (sTrig.IsTimeOverThresholdDeconvoluted() && !fSilenceTOTdStations) ||
            (sTrig.IsMultiplicityOfPositiveSteps() && !fSilenceMoPSStations))
          f3TOTMixCandidateStations.push_back(sp);
      }
    }

    const bool is3TOT = IsSEvent3TOT(f3TOTCandidateStations);
    const bool is3TOTd = IsSEvent3TOT(f3TOTdCandidateStations);
    const bool is3MoPS = IsSEvent3TOT(f3MoPSCandidateStations);
    const bool is3TOTMix = IsSEvent3TOT(f3TOTMixCandidateStations); // superset of is3TOT is3TOTd is3MoPS

    if (is3TOT)
      ++namedCounters["SdEventSelector/3TOT"];
    if (is3TOTd)
      ++namedCounters["SdEventSelector/3TOTd"];
    if (is3MoPS)
      ++namedCounters["SdEventSelector/3MoPS"];
    if (is3TOTMix)
      ++namedCounters["SdEventSelector/3TOTMix"];

    if (fEnableBottomUpSelection && !BottomUpSelection(is3TOTMix, event)) {
      INFO("Bottom-up-selection failed!");
      ++namedCounters["SdEventSelector/BottomUpSelectionFailed"];
      return eContinueLoop;
    }

    {
      int lonely = 0;
      while (RejectLonelyStations())
        ++lonely;
      if (lonely)
        ++namedCounters["SdEventSelector/LonelyStationEvents"];
    }

    ostringstream info;

    int t4TriggerPattern = 0;
    if (fEnableT4Trigger) {

      const bool is4C1 = IsSEvent4C1();
      const bool isT4 = is3TOTMix || is4C1;

      info << "T4 = "
           << GetT4TriggerName(isT4, isFD, is3TOT, is3TOTd, is3MoPS, is3TOTMix, is4C1);

      if (is4C1)
        ++namedCounters["SdEventSelector/4C1"];
      if (isT4)
        ++namedCounters["SdEventSelector/T4"];

      if (fRejectNonT4Events && !isT4) {
        INFO(info);
        return eContinueLoop;
      }

      t4TriggerPattern =
        (isFD      * ShowerSRecData::eT4_FD)      |
        (is3TOT    * ShowerSRecData::eT4_3TOT)    |
        (is3TOTd   * ShowerSRecData::eT4_3TOTd)   |
        (is3MoPS   * ShowerSRecData::eT4_3MoPS)   |
        (is3TOTMix * ShowerSRecData::eT4_3TOTMix) |
        (is4C1     * ShowerSRecData::eT4_4C1);

      shSRec.SetT4Trigger(t4TriggerPattern);
      shSRec.SetReconstructionSeed(fSeedStations);

    }

    if (fEnableT5Trigger) {

      fCandidateStations.clear();

      for (auto& s : sEvent.CandidateStationsRange()) {
        const auto& ds = sDet.GetStation(s);
        fCandidateStations.emplace_back(&s, &ds);
      }

      const auto t5Stations = CalculateT5Trigger(sEvent);
      const unsigned int hottestStation = t5Stations.first;
      const int nActive = t5Stations.second.size();
      const bool t5 = (hottestStation && nActive >= fMinNumberOfActiveStations);

      // Fill vector denoting which T5 neighbor stations were UUBs
      vector<int> t5NeighborsUUB;
      for (const auto sid : t5Stations.second) {
        const bool isUUB = sDet.GetStation(sid).IsUUB();
        if (isUUB)
          t5NeighborsUUB.push_back(sid);
      }

      info << (info.str().length() ? "; " : "")
           << "T5 = " << t5 << ", "
           << nActive << "T5, "
              "stations = " << hottestStation << "(hottest)";
      for (const auto id : t5Stations.second) {
        const auto& s = sEvent.GetStation(id);
        info << ' ' << id << '('
             << (s.IsCandidate() ?
                   "candidate" :
                   (s.IsSilent() ?
                      "silent" :
                      ("rejected:" +
                       AsBinary(s.GetRejectionStatus(), StationConstants::eNumRejectionStatusBits))
                   )
                ) << ')';
      }

      if (t5)
        ++namedCounters["SdEventSelector/T5"];

      if (fRejectNonT5Events && !t5) {
        INFO(info);
        return eContinueLoop;
      }

      shSRec.SetT5Trigger(t5*ShowerSRecData::eT5_Prior);
      shSRec.SetT5HottestStation(hottestStation);
      shSRec.SetT5PriorActiveNeighbors(t5Stations.second);
      shSRec.SetT5PriorActiveUUBNeighbors(t5NeighborsUUB);

    }

    const int id = GetBadPeriod(eventTime);
    if (id) {
      event.GetRecShower().GetSRecShower().SetBadPeriodId(id);
      ++namedCounters["SdEventSelector/BadPeriodEvents"];
      info << "; Bad period id=" << id;
      INFO(info);
      if (fRejectBadPeriods)
        return eContinueLoop;
    }

    if (!info.str().empty())
      INFO(info);

    ++namedCounters["SdEventSelector/SelectedEvent"];

    return eSuccess;
  }


  VModule::ResultFlag
  SdEventSelector::Finish()
  {
    return eSuccess;
  }


  int
  SdEventSelector::GetBadPeriod(const TimeStamp& time)
    const
  {
    if (fBadPeriods.empty())
      return 0;

    // gives first period not less than time
    auto low =
      lower_bound(fBadPeriods.begin(), fBadPeriods.end(), TimeRange(time, time));

    if (low == fBadPeriods.begin())
      return 0;

    // check if the last period less than time contains time
    --low;
    return low->IsInRange(time) ? low->GetId() : 0;
  }


  /*
    to find 3TOT we need a triangle with (for shells see the table at top):
    - 3 edges from shell 1 (basic equilateral triangle) or
    - 2 edges from shell 1 and one from shell 2 (stretched triangle)
  */
  bool
  SdEventSelector::IsSEvent3TOT(const StationPairList& candidateStations)
    const
  {
    if (candidateStations.size() < 3)
      return false;

    const auto spBegin = candidateStations.begin();
    for (auto spIt = spBegin, spEnd = candidateStations.end(); spIt != spEnd; ++spIt) {

      const auto& si = *spIt->first;
      const auto& dsi = *spIt->second;

      for (auto spJt = spBegin; spJt != spEnd; ++spJt) {
        if (spJt == spIt)
          continue;

        const auto& sj = *spJt->first;
        const auto& dsj = *spJt->second;

        // require i-j ~ 1500 m
        const double distij = Distance(dsi, dsj);
        if (distij < fFirstCrownDistanceRange[0] || distij > fFirstCrownDistanceRange[1] ||
            distij < CTimeDifferenceMargin(si, sj))
          continue;

        for (auto spKt = spBegin; spKt != spEnd; ++spKt) {
          if (spKt == spJt || spKt == spIt)
            continue;

          const auto& sk = *spKt->first;
          const auto& dsk = *spKt->second;

          // require i-k ~ 1500 m
          const double distik = Distance(dsi, dsk);
          if (distik < fFirstCrownDistanceRange[0] || distik > fFirstCrownDistanceRange[1] ||
              distik < CTimeDifferenceMargin(si, sk))
            continue;

          // require j-k ~ 1500 or 2600 m
          const double distjk = Distance(dsj, dsk);
          if (distjk < fFirstCrownDistanceRange[0] || distjk > fSkewMaximumDistance ||
              distjk < CTimeDifferenceMargin(sj, sk))
            continue;

          // Found a 3TOT
          return true;
        }
      }
    }

    return false;
  }


  /*
    for 4C1 configuration we have to find a station (si) with 3 neighbors from
    shell 1 (sj, sk, sl)
  */
  bool
  SdEventSelector::IsSEvent4C1()
    const
  {
    if (f4C1CandidateStations.size() < 4)
      return false;

    const auto spBegin = f4C1CandidateStations.begin();
    for(auto spIt = spBegin, spEnd = f4C1CandidateStations.end(); spIt != spEnd; ++spIt) {

      const auto& si = *spIt->first;
      const auto& dsi = *spIt->second;

      for (auto spJt = spBegin; spJt != spEnd; ++spJt) {
        if (spJt == spIt)
          continue;

        const auto& sj = *spJt->first;
        const auto& dsj = *spJt->second;

        // require i-j ~ 1500 m
        const double distij = Distance(dsi, dsj);
        if (distij < fFirstCrownDistanceRange[0] || distij > fFirstCrownDistanceRange[1] ||
            distij < CTimeDifferenceMargin(si, sj))
          continue;

        for (auto spKt = Next(spJt); spKt != spEnd; ++spKt) {
          if (spKt == spIt)
            continue;

          const auto& sk = *spKt->first;
          const auto& dsk = *spKt->second;

          // require i-k ~ 1500 m
          const double distik = Distance(dsi, dsk);
          if (distik < fFirstCrownDistanceRange[0] || distik > fFirstCrownDistanceRange[1] ||
              distik < CTimeDifferenceMargin(si, sk))
            continue;

          if (Distance(dsj, dsk) < CTimeDifferenceMargin(sj, sk))
            continue;

          for (auto spLt = Next(spKt); spLt != spEnd; ++spLt) {
            if (spLt == spIt)
              continue;

            const auto& sl = *spLt->first;
            const auto& dsl =*spLt->second;

            // require i-l ~ 1500 m
            const double distil = Distance(dsi, dsl);
            if (distil < fFirstCrownDistanceRange[0] || distil > fFirstCrownDistanceRange[1] ||
                distil < CTimeDifferenceMargin(si, sl))
              continue;

            if (Distance(dsj, dsl) < CTimeDifferenceMargin(sj, sl))
              continue;

            if (Distance(dsk, dsl) < CTimeDifferenceMargin(sk, sl))
              continue;

            // found 4C1
            return true;
          }
        }
      }
    }

    return false;
  }


  pair<int, vector<int>> // hottest station id, number of active neighbors
  SdEventSelector::CalculateT5Trigger(const SEvent& sEvent)
    const
  {
    const Station* maxStation = nullptr;
    {
      double maxSignal = 0;
      for (const auto& s : sEvent.CandidateStationsRange()) {
        const double signal = s.GetRecData().GetTotalSignal();
        if (maxSignal < signal) {
          maxSignal = signal;
          maxStation = &s;
        }
      }

      if (!maxStation) {
        INFO("No station with maximal signal found");
        return make_pair(0U, vector<int>());
      }
    }

    const auto& sDet = det::Detector::GetInstance().GetSDetector();
    const auto maxPos = sDet.GetStation(*maxStation).GetPosition();

    const auto badBits = ~StationConstants::kGoodBitsT5;

    pair<unsigned int, vector<int>> ret;
    ret.first = maxStation->GetId();
    auto& functioningNeighbors = ret.second;
    unsigned int nCandidates = 0;
    unsigned int nSilent = 0;
    unsigned int nRejected = 0;
    for (const auto& s : sEvent.StationsRange()) {
      if (!s.IsRejected() || !(s.GetRejectionStatus() & badBits)) {
        const auto& pos = sDet.GetStation(s).GetPosition();
        const double dist = Distance(maxPos, pos);
        if (fFirstCrownDistanceRange[0] < dist && dist < fFirstCrownDistanceRange[1]) {
          functioningNeighbors.push_back(s.GetId());
          nCandidates += s.IsCandidate();
          nSilent += s.IsSilent();
          nRejected += s.IsRejected();
        }
      }
    }

    if (int(functioningNeighbors.size()) < fMinNumberOfActiveStations) {
      ostringstream info;
      info << "Only " << functioningNeighbors.size() << " "
              "(out of " << fMinNumberOfActiveStations << ") "
              "functioning neighbors for max-signal station "
           << maxStation->GetId() << " found "
              "(candidates: " << nCandidates << ", "
               "silent: " << nSilent << ", "
               "rejected: " << nRejected << ").";
      INFO(info);
    }

    return ret;
  }


  bool
  SdEventSelector::BottomUpSelection(const bool is3TOTMix, Event& event)
  {
    // must have valid f3TOTMixCandidateStations!

    set<StationPair> all4C1Stations;

    if (!is3TOTMix && f4C1CandidateStations.size() >= 4) {
      // repeat 4C1 search and find ALL stations that satisfy 4C1 condition

      const auto spBegin = f4C1CandidateStations.begin();
      for (auto spIt = spBegin, spEnd = f4C1CandidateStations.end(); spIt != spEnd; ++spIt) {

        const auto& si = *spIt->first;
        const auto& dsi = *spIt->second;

        for (auto spJt = spBegin; spJt != spEnd; ++spJt) {
          if (spJt == spIt)
            continue;

          const auto& sj = *spJt->first;
          const auto& dsj = *spJt->second;

          // require i-j ~ 1500 m
          const double distij = Distance(dsi, dsj);
          if (distij < fFirstCrownDistanceRange[0] || distij > fFirstCrownDistanceRange[1] ||
              distij < CTimeDifferenceMargin(si, sj))
            continue;

          for (auto spKt = Next(spJt); spKt != spEnd; ++spKt) {
            if (spKt == spIt)
              continue;

            const auto& sk = *spKt->first;
            const auto& dsk = *spKt->second;

            // require i-k ~ 1500 m
            const double distik = Distance(dsi, dsk);
            if (distik < fFirstCrownDistanceRange[0] || distik > fFirstCrownDistanceRange[1] ||
                distik < CTimeDifferenceMargin(si, sk))
              continue;

            if (Distance(dsj, dsk) < CTimeDifferenceMargin(sj, sk))
              continue;

            for (auto spLt = Next(spKt); spLt != spEnd; ++spLt) {
              if (spLt == spIt)
                continue;

              const auto& sl = *spLt->first;
              const auto& dsl =*spLt->second;

              // require i-l ~ 1500 m
              const double distil = Distance(dsi, dsl);
              if (distil < fFirstCrownDistanceRange[0] || distil > fFirstCrownDistanceRange[1] ||
                  distil < CTimeDifferenceMargin(si, sl))
                continue;

              if (Distance(dsj, dsl) < CTimeDifferenceMargin(sj, sl))
                continue;

              if (Distance(dsk, dsl) < CTimeDifferenceMargin(sk, sl))
                continue;

              all4C1Stations.insert(*spIt);
              all4C1Stations.insert(*spJt);
              all4C1Stations.insert(*spKt);
              all4C1Stations.insert(*spLt);
            } // for spLt
          } // for spKt
        } // for spJt
      } // for spIt

    } // if not 3TOT

    fSeedStations.clear();

    // add stations above
    for (const auto& sp : all4C1Stations) {
      const auto e = f3TOTMixCandidateStations.end();
      if (find(f3TOTMixCandidateStations.begin(), e, sp) == e)
        f3TOTMixCandidateStations.push_back(sp);
    }

    if (f3TOTMixCandidateStations.size() < 3) {
      INFO("No seed found.");
      return false;
    }

    double maxSignalSum = 0;
    StationPair seed[3] = { StationPair(0, 0) };

    // find triangle (3TOT && 4C1 above) seed stations
    const auto spBegin = f3TOTMixCandidateStations.begin();
    for (auto spIt = spBegin, spEnd = f3TOTMixCandidateStations.end(); spIt != spEnd; ++spIt) {

      const auto& si = *spIt->first;
      const auto& dsi = *spIt->second;

      for (auto spJt = spBegin; spJt != spEnd; ++spJt) {

        if (spJt == spIt)
          continue;

        const auto& sj = *spJt->first;
        const auto& dsj = *spJt->second;

        const double distij = Distance(dsi, dsj);

        // require i-j ~ 1500 m
        if (distij < fFirstCrownDistanceRange[0] || distij > fFirstCrownDistanceRange[1] ||
            distij < CTimeDifferenceMargin(si, sj))
          continue;

        for (auto spKt = spBegin; spKt != spEnd; ++spKt) {

          if (spKt == spIt || spKt == spJt)
            continue;

          const auto& sk = *spKt->first;
          const auto& dsk = *spKt->second;

          const double distik = Distance(dsi, dsk);

          // require i-k ~ 1500 m
          if (distik < fFirstCrownDistanceRange[0] || distik > fFirstCrownDistanceRange[1] ||
              distik < CTimeDifferenceMargin(si, sk))
            continue;

          const double distjk = Distance(dsj, dsk);

          if (distjk < fFirstCrownDistanceRange[0] || distjk > fSkewMaximumDistance ||
              distjk < CTimeDifferenceMargin(sj, sk))
            continue;

          const double sigi = si.GetRecData().GetTotalSignal();
          const double sigj = sj.GetRecData().GetTotalSignal();
          const double sigk = sk.GetRecData().GetTotalSignal();
          const double signalSum = sigi + sigj + sigk;

          if (signalSum > maxSignalSum) {
            maxSignalSum = signalSum;
            // find max
            if (sigi > sigj) {
              if (sigi > sigk) {
                seed[0] = *spIt;
                seed[1] = *spJt;
                seed[2] = *spKt;
              } else {
                seed[0] = *spKt;
                seed[1] = *spIt;
                seed[2] = *spJt;
              }
            } else {
              if (sigj > sigk) {
                seed[0] = *spJt;
                seed[1] = *spIt;
                seed[2] = *spKt;
              } else {
                seed[0] = *spKt;
                seed[1] = *spIt;
                seed[2] = *spJt;
              }
            }
          }

        } // for spKt
      } // for spJt
    } // for spIt

    if (!seed[0].first) {

      INFO("No seed found.");
      return false;

    } else {

      fSeedStations.resize(3);
      for (int i = 0; i < 3; ++i)
        fSeedStations[i] = seed[i].first->GetId();

      const auto& s1 = *seed[0].first;
      const auto& t1 = s1.GetRecData().GetSignalStartTime();
      const auto& x1 = seed[0].second->GetPosition();

      const auto& s2 = *seed[1].first;
      const double ct12 =
        kSpeedOfLight*(t1 - s2.GetRecData().GetSignalStartTime()).GetInterval();
      const auto x21 = seed[1].second->GetPosition() - x1;

      const auto& s3 = *seed[2].first;
      const double ct13 =
        kSpeedOfLight*(t1 - s3.GetRecData().GetSignalStartTime()).GetInterval();
      const auto x31 = seed[2].second->GetPosition() - x1;

      const double sx21 = x21.GetMag();
      const auto vi = x21 / sx21;
      const double taui = ct12/sx21;
      const double sx31 = x31.GetMag();
      const auto vj = x31 / sx31;
      const double tauj = ct13/sx31;
      const double ij = vi*vj;
      const double d = 1 - ij*ij;
      const double gamma2 = (1 - (taui*taui - 2*ij*taui*tauj + tauj*tauj)/d)/d;

      const auto& siteCS = det::Detector::GetInstance().GetSiteCoordinateSystem();

      Vector axis;

      if (gamma2 >= 0) {

        const double alpha = (taui - tauj*ij) / d;
        const double beta = (tauj - taui*ij) / d;
        const double gamma = sqrt(gamma2);
        Vector vk = Cross(vi, vj);
        // flip to zenith side
        if (vk.GetZ(siteCS) < 0)
          vk *= -1;
        axis = alpha*vi + beta*vj + gamma*vk;

      } else {

        // fit
        INFO("Unphysical seed, trying minimization.");

        TMinuit minuit(1);
        int errFlag = 0;
        double argList[10];
        argList[0] = -1;
        minuit.mnexcm("SET PRINTOUT", argList, 1, errFlag);
        minuit.mnexcm("SET NOWARNINGS", argList, 0, errFlag);
        minuit.SetPrintLevel(-1);

        const auto vk = Normalized(vj - ij * vi);

        minuit.SetFCN(SeedFit::Chi2);
        SeedFit::gScalarProductIJ = ij;
        const double jk = vj * vk;
        SeedFit::gScalarProductJK = jk;
        SeedFit::gTimeI = taui;
        SeedFit::gTimeJ = tauj;

        argList[0] = 1;
        minuit.mnexcm("SET ERRORDEF", argList, 1, errFlag);

        const double phi0 = atan2(tauj*jk, taui + tauj*ij);
        minuit.mnparm(0, "phi", phi0, 0.5, 0, 0, errFlag);

        argList[0] = 500;
        minuit.mnexcm("MINIMIZE", argList, 1, errFlag);
        if (errFlag) {
          minuit.mnexcm("MINOS", argList, 0, errFlag);
          if (errFlag)
            return false;
        }

        double phi;
        double phiErr;
        minuit.GetParameter(0, phi, phiErr);

        const double u = cos(phi);
        const double v = sin(phi);
        axis = u * vi + v * vk;
      }

      {
        ostringstream info;
        info << "Seed theta = " << axis.GetTheta(siteCS)/degree << " deg, "
                "phi = " << axis.GetPhi(siteCS)/degree << " deg (site)";
        INFO(info);
      }

      auto& ssRec = GetSRecShower(event);
      ssRec.SetSeedAxis(axis);
      ssRec.SetSeedBarycenter(x1 + (1/3.)*(x21 + x31));
      // find accidentals
      map<int, double> rejected;
      for (auto spIt = fCandidateStations.begin(); spIt != fCandidateStations.end(); ) {
        auto& s = *spIt->first;
        const double delay =
          (s.GetRecData().GetSignalStartTime() - t1).GetInterval() +
          (axis*(spIt->second->GetPosition() - x1)) / kSpeedOfLight;

        s.SetBottomUpResidual(delay);

        // is delay within [Early, Late]
        const auto& sTrig = s.GetTriggerData();
        if ((sTrig.IsTimeOverThreshold() && IsInRange(delay, fSeedPlaneFrontTimeResidualRangeToT)) ||
            (sTrig.IsTimeOverThresholdDeconvoluted() && IsInRange(delay, fSeedPlaneFrontTimeResidualRangeToTd)) ||
            (sTrig.IsMultiplicityOfPositiveSteps() && IsInRange(delay, fSeedPlaneFrontTimeResidualRangeMoPS)) ||
            (sTrig.IsT2Threshold() && IsInRange(delay, fSeedPlaneFrontTimeResidualRangeTh2)) ||
            (sTrig.IsT1Threshold() && IsInRange(delay, fSeedPlaneFrontTimeResidualRangeTh1)))
          ++spIt;
        else {
          s.SetRejected(StationConstants::eOutOfTime);
          rejected[s.GetId()] = delay;
          f3TOTMixCandidateStations.remove(*spIt);
          f4C1CandidateStations.remove(*spIt);
          spIt = fCandidateStations.erase(spIt);
        }
      }

      if (!rejected.empty()) {
        ostringstream info;
        info << "Station" << String::Plural(rejected)
             << " not compatible with seed reconstruction:\n";
        TabularStream tab("r|.");
        tab <<              endc << "time" << endr
            << "station" << endc << "diff [ns]" << endr << hline;
        for (const auto& r : rejected)
          tab << r.first << ' ' << endc << ' ' << r.second/nanosecond << endr;
        tab << delr;
        info << tab;
        INFO(info);
      }

    }

    return true;
  }


  bool
  SdEventSelector::RejectLonelyStations()
  {
    vector<int> ssLonely;

    for (auto spIt = fCandidateStations.begin(); spIt != fCandidateStations.end(); ) {
      const auto& dsi = *spIt->second;
      int closeNeighbors = 0;
      int remoteNeighbors = 0;
      for (auto spJt = fCandidateStations.begin(); spJt != fCandidateStations.end(); ++spJt) {
        if (spJt == spIt)
          continue;
        const double distij = Distance(dsi, *spJt->second);
        if (distij < fLonelyIfNoneInDistance)
          ++closeNeighbors;
        if (distij < fLonelyIfOneInDistance)
          ++remoteNeighbors;
      }
      if (!closeNeighbors || remoteNeighbors <= 1) {
        auto& s = *spIt->first;
        ssLonely.push_back(s.GetId());
        s.SetRejected(StationConstants::eLonely);
        f3TOTMixCandidateStations.remove(*spIt);
        f4C1CandidateStations.remove(*spIt);
        spIt = fCandidateStations.erase(spIt);
      } else
        ++spIt;
    }

    if (!ssLonely.empty()) {
      ostringstream info;
      info << String::OfSortedIds(ssLonely) << " rejected as lonely.";
      INFO(info);
      return true;
    }

    return false;
  }

}