SdEventPosteriorSelector.cc

Go to the documentation of this file.

#include "SdEventPosteriorSelector.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/SEventConstants.h>
#include <sevt/Header.h>
#include <sevt/EventTrigger.h>
#include <sevt/Station.h>
#include <sevt/StationTriggerData.h>
#include <sevt/StationRecData.h>
#include <sevt/SortCriteria.h>

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

#include <boost/range/iterator_range.hpp>

using namespace fwk;
using namespace evt;
using namespace sevt;
using namespace utl;
using namespace std;


namespace SdEventPosteriorSelectorOG {

  inline
  bool
  IsAlongside(const Vector& p1, const Vector& p2, const Vector& p3)
  {
    return Cross(p1, p2) * Cross(p1, p3) >= 0;
  }


  // projects core along axis to a plane formed by vectors b-a and c-a given triangle point a, b, c
  inline
  Point
  ProjectPointAlongAxisToTriangle(const Point& core, const Vector& axis,
                                  const Point& a, const Point& b, const Point& c)
  {
    const Vector n = Cross(b - a, c - a);
    return core + ((n*(a - core)) / (n*axis)) * axis;
  }


  inline
  bool
  IsInsideTriangle(const Point& p, const Point& a, const Point& b, const Point& c)
  {
    return IsAlongside(a-b, p-b, c-b) &&
           IsAlongside(b-c, p-c, a-c) &&
           IsAlongside(c-a, p-a, b-a);
  }


  pair<vector<int>, unsigned int>
  SdEventPosteriorSelector::CalculateT5PosteriorTrigger(const Event& event)
    const
  {
    const auto& sEvent = event.GetSEvent();
    const auto& shRec = event.GetRecShower();
    if (!shRec.HasSRecShower() || sEvent.CandidateStationsBegin() == sEvent.CandidateStationsEnd())
      return make_pair(vector<int>(), SEventConstants::eNone);
    const auto& shSRec = shRec.GetSRecShower();

    const Station* hottestStation = nullptr;
    if (const unsigned int id = shSRec.GetT5HottestStation())
      hottestStation = &sEvent.GetStation(id);
    else {
      sEvent.SortStations(ByDecreasingSignal());
      hottestStation = &*sEvent.CandidateStationsBegin();
    }

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

    const auto& core = shSRec.GetCorePosition();

    const auto badBits = ~StationConstants::kGoodBitsT5;

    // find crown C1 around the hottest station, and stations for triangle calculation below
    vector<Point> neighborPos;
    pair<vector<int>, unsigned int> result;
    result.second = SEventConstants::eNone;
    auto& neighbors = result.first;
    for (const auto& s : sEvent.StationsRange()) {
      if (s.IsRejected() && (s.GetRejectionStatus() & badBits)) // make sure no bad bit is set
        continue;
      const auto& pos = sDet.GetStation(s).GetPosition();
      const double dist = Distance(hotPos, pos);
      if (fFirstCrownDistanceRange[0] < dist && dist < fFirstCrownDistanceRange[1]) {
        neighbors.push_back(s.GetId());
        neighborPos.push_back(pos);
      }
    }

    if (neighbors.size() > 6) {
      const string msg = "Number of C1 stations should not be larger than 6!";
      if (fRejectNonT5PosteriorEvents) {
        ERROR(msg);
        throw OutOfBoundException(msg);
      } else
        WARNING(msg);
    }

    if (neighborPos.size() < 2)
      return result;

    const auto& axis = shSRec.GetAxis();

    unsigned int nTriangles[2] = { 0, 0 };
    bool isCoreClosestToHottest = false;
    bool isCoreInHottestVoronoi = false;
    const int n = neighborPos.size();
    for (int i = 0, n1 = n - 1; i < n1; ++i) {
      for (int j = i+1; j < n; ++j) {
        const auto& ci = neighborPos[i];
        const auto& cj = neighborPos[j];
        const double dij = Distance(ci, cj);
        // consider only nearest and next-nearest neighbors in C1
        if (dij > fSkewMaximumDistance)
          continue;
        const auto pc = ProjectPointAlongAxisToTriangle(core, axis, hotPos, ci, cj);
        if (!IsInsideTriangle(pc, hotPos, ci, cj))
          continue;
        const bool isIsosceles = (dij > fFirstCrownDistanceRange[1]);
        ++nTriangles[isIsosceles];
        {
          const double dh = (pc - hotPos).GetMag2();
          if (dh < (pc - ci).GetMag2() && dh < (pc - cj).GetMag2())
            isCoreClosestToHottest = true;
        }
        {
          const Line shower(pc, axis);
          const double dh = Distance(shower, hotPos);
          if (dh < Distance(shower, ci) && dh < Distance(shower, cj))
            isCoreInHottestVoronoi = true;
        }
      }
    }
    result.second =
      (bool(nTriangles[0]) * SEventConstants::eEquilateral) |
      (bool(nTriangles[1]) * SEventConstants::eIsosceles) |
      (isCoreClosestToHottest * SEventConstants::eCoreClosestToHottest) |
      (isCoreInHottestVoronoi * SEventConstants::eCoreInHottestVoronoi);
    return result;
  }


  VModule::ResultFlag
  SdEventPosteriorSelector::Init()
  {
    Branch topB = CentralConfig::GetInstance()->GetTopBranch("SdEventSelector");
    topB.GetChild("EnableT5PosteriorTrigger").GetData(fEnableT5Trigger);
    topB.GetChild("RejectNonT5PosteriorEvents").GetData(fRejectNonT5PosteriorEvents);
    topB.GetChild("MinNumberOfActiveStationsPosterior").GetData(fMinNumberOfActiveStationsPosterior);
    topB.GetChild("FirstCrownDistanceRange").GetData(fFirstCrownDistanceRange);
    topB.GetChild("SkewMaximumDistance").GetData(fSkewMaximumDistance);
    return eSuccess;
  }


  VModule::ResultFlag
  SdEventPosteriorSelector::Run(evt::Event& event)
  {
    if (!event.HasSEvent())
      return eContinueLoop;
    const auto& sEvent = event.GetSEvent();

    if (!event.HasRecShower())
      return eContinueLoop;
    auto& shRec = event.GetRecShower();
    
    if (!shRec.HasSRecShower())
      return eContinueLoop;
    auto& shSRec = shRec.GetSRecShower();

    if (fEnableT5Trigger) {
      const auto t5NeighborsTriangle = CalculateT5PosteriorTrigger(event);
      const int nActive = t5NeighborsTriangle.first.size();
      const unsigned int triangleType = t5NeighborsTriangle.second;
      const bool t5Posterior = (nActive >= fMinNumberOfActiveStationsPosterior && triangleType);

      ostringstream info;
      info << "T5Posterior = " << (t5Posterior ? "true" : "false") << ", "
              "NT5 = "  << nActive << ", "
              "type = " << triangleType << ':'
           << (triangleType ? (
                string((triangleType & SEventConstants::eCoreInHottestVoronoi) ? " voronoi" : "") +
                string((triangleType & SEventConstants::eCoreClosestToHottest) ? " closest" : "") +
                string((triangleType & SEventConstants::eIsosceles) ? " isosceles" : "") +
                string((triangleType & SEventConstants::eEquilateral) ? " equilateral" : "")
              ) : " none") << ", "
              "station ids = " << shSRec.GetT5HottestStation() << "(hottest)";
      for (const auto id : t5NeighborsTriangle.first) {
        const auto& s = sEvent.GetStation(id);
        info << ' ' << id << '('
             << (s.IsCandidate() ?
                   "candidate" :
                   (s.IsSilent() ?
                      "silent" :
                      ("rejected:" + AsBinary(s.GetRejectionStatus(), StationConstants::eNumRejectionStatusBits))
                   )
                ) << ')';
      }
      INFO(info);

      ++RunController::GetInstance().GetRunData().GetNamedCounters()["SdEventPosteriorSelector/T5Post"];

      if (fRejectNonT5PosteriorEvents && !t5Posterior)
        return eContinueLoop;

      int t5Trigger = shSRec.GetT5Trigger();
      info.str("");
      info << "t5trigger " << t5Trigger << " -> ";
      AsBitArray(t5Trigger).Mask(ShowerSRecData::eT5_Posterior, t5Posterior);
      shSRec.SetT5Trigger(t5Trigger);
      info << t5Trigger;
      INFO(info);

      shSRec.SetT5PosteriorActiveNeighbors(t5NeighborsTriangle.first);
      shSRec.SetT5PosteriorCoreTriangle(triangleType);
    }

    return eSuccess;
  }

}