RdClusterFinder.cc

Go to the documentation of this file.

#include "RdClusterFinder.h"

#include <fwk/CentralConfig.h>
#include <utl/Branch.h>
#include <utl/RadioGeometryUtilities.h>

#include <det/Detector.h>
#include <rdet/RDetector.h>

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

#include <revt/REvent.h>
#include <revt/StationRecData.h>

#include <vector>
#include <algorithm>

using namespace std;
using namespace fwk;
using namespace utl;
using namespace evt;
using namespace revt;


namespace RdClusterFinder {

  bool
  SortPairByDistance(const pair<double, int> x, const pair<double, int> y)
  {
    return (x.first > y.first);
  }


  VModule::ResultFlag
  RdClusterFinder::Init()
  {
    Branch topBranch = CentralConfig::GetInstance()->GetTopBranch("RdClusterFinder");
    topBranch.GetChild("InfoLevel").GetData(fInfoLevel);
    topBranch.GetChild("MaximumDiscontinuity").GetData(fMaximumDiscontinuity);
    topBranch.GetChild("RejectLonelyStations").GetData(fRejectLonelyStations);
    topBranch.GetChild("LonelyIfNoneInDistance").GetData(fLonelyIfNoneInDistance);
    topBranch.GetChild("LonelyIfOneInDistance").GetData(fLonelyIfOneInDistance);
    topBranch.GetChild("LonelyIfNoneInDistancePhase3").GetData(fLonelyIfNoneInDistancePhase3);
    topBranch.GetChild("LonelyIfOneInDistancePhase3").GetData(fLonelyIfOneInDistancePhase3);
    topBranch.GetChild("IgnoreRejectedStations").GetData(fIgnoreRejectedStations);
    topBranch.GetChild("IgnoreExcludedStations").GetData(fIgnoreExcludedStations);

    if (fIgnoreRejectedStations > 2)
      ERROR("Wrong config. IgnoreRejectedStations must be 0, 1 or 2.");

    return eSuccess;
  }


  VModule::ResultFlag
  RdClusterFinder::Run(Event& event)
  {
    // Check if there are no events at all
    if (!event.HasREvent()) {
      WARNING("eContinueLoop, No radio event found!");
      return eContinueLoop;
    }

    REvent& rEvent = event.GetREvent();
    ostringstream info;

    // cluster finder will not work properly for less than 3 signal stations and reject all remaining stations
    if (rEvent.GetNumberOfSignalStations() < 3) {
      return eSuccess;
    }

    const rdet::RDetector& rDetector = det::Detector::GetInstance().GetRDetector();
    vector<int> selectedCluster;
    vector<int> deselectedStations;

    const Point& showerCore = event.GetRecShower().GetRRecShower().GetReferenceCorePosition(event);
    const utl::Vector& showerAxis = event.GetRecShower().GetRRecShower().GetReferenceAxis(event);

    // distance between Sd axis and signal stations and station id
    vector<pair<double, int> > distanceToSdAxis;
    for (const auto& station : rEvent.StationsRange()) {

      if (!UseStation(station))
        continue;

      const Point& stationPosition = rDetector.GetStation(station).GetPosition();
      const float distance = utl::RadioGeometryUtilities::GetDistanceToAxis(showerAxis, showerCore,
                                                                            stationPosition);
      distanceToSdAxis.emplace_back(distance, station.GetId());
    }
    sort(distanceToSdAxis.begin(), distanceToSdAxis.end(), SortPairByDistance);

    // check axis distances for jumps larger than fMaximumDiscontinuity and deselect
    // loop from furthest away from axis to 2nd closest stations
    for (auto itCurrent = distanceToSdAxis.begin(); itCurrent != prev(distanceToSdAxis.end()); ++itCurrent) {
      selectedCluster.push_back(itCurrent->second);

      auto itNext = next(itCurrent);
      const double deltaR = itCurrent->first - itNext->first;
      if (deltaR > fMaximumDiscontinuity) {
        copy(selectedCluster.begin(), selectedCluster.end(), std::back_inserter(deselectedStations));
        selectedCluster.clear();
      }
    }

    // reject the deselected stations
    for (const auto& deselectedStation : deselectedStations) {
      Station& theStation = rEvent.GetStation(deselectedStation);
      theStation.SetRejectedReason(revt::eLonelyStation);

      //Display some information for low level bug tracking
      info.str("");
      info << "Deselected station ID : " << deselectedStation;
      INFOFinal(info.str());
    }

    // get reference coordinate system of detector (usually PampaAmarilla)
    // needed to detect AERA phase 3 via position aka south
    const CoordinateSystemPtr referenceCS =
      det::Detector::GetInstance().GetReferenceCoordinateSystem();

    /* lonely station rejection
       count close and remote number of neighboring stations that are interpreted as signal station.
       iStation: current station under investigation, only iterate signal stations
       jStation: all possible neighbour stations, use all station iterator to include rejected
                 and excluded stations with may be counted as a signal station depending on xml config */
    if (fRejectLonelyStations) {
      for (auto& iStation : rEvent.SignalStationsRange()) {
        unsigned int closeNeighbors = 0;
        unsigned int remoteNeighbors = 0;

        double closeDist = 0;
        double remoteDist = 0;

        /* check for iStation in phase 3 and use larger distances
         We could hardcode a list of phase 3 IDs but IDs can change if station is moved.
         Everything south of 14.6 km im PampaAmarilla is considered phase 3 */
        if (rDetector.GetStation(iStation.GetId()).GetPosition().GetY(referenceCS) <  14.6 * km) {
          closeDist = fLonelyIfNoneInDistancePhase3;
          remoteDist = fLonelyIfOneInDistancePhase3;
        } else {
          closeDist = fLonelyIfNoneInDistance;
          remoteDist = fLonelyIfOneInDistance;
        }

        for (const auto& jStation : rEvent.AllStationsRange()) {

          if (iStation.GetId() == jStation.GetId())
            continue;

          if (!UseStation(jStation))
            continue;

          info.str("");
          info << "checking stations: " << iStation.GetId() << ", " << jStation.GetId();
          INFODebug(info.str());

          const double distij = (rDetector.GetStation(iStation.GetId()).GetPosition() -
                                 rDetector.GetStation(jStation.GetId()).GetPosition()).GetMag();

          if (distij < closeDist)
            ++closeNeighbors;
          if (distij < remoteDist)
            ++remoteNeighbors;
        }

        if (!closeNeighbors || remoteNeighbors <= 1) {
          iStation.SetRejectedReason(revt::eLonelyStation);

          info.str("");
          if (!closeNeighbors)
            info << "Station with ID " << iStation.GetId()
                 << " rejected as lonely due to # of close neighbors = " << closeNeighbors;
          else if (remoteNeighbors <= 1)
            info << "Station with ID " << iStation.GetId()
                 << " rejected as lonely due to # of remote neighbors = " << remoteNeighbors;
          INFOFinal(info.str());
        }
      }
    }

    return eSuccess;
  }


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


  bool
  RdClusterFinder::UseStation(const revt::Station& station)
  {
    if (station.IsExcluded() && !fIgnoreExcludedStations)
      return false;

    if (!CheckRejection(station))
      return false;

    // if fIgnoreExcludedStations they are skipped already above, if not they count as signal
    if (!(station.IsExcluded() || station.GetRecData().GetPulseFound()))
      return false;

    return true;
  }


  bool
  RdClusterFinder::CheckRejection(const revt::Station& station)
  {
    switch (fIgnoreRejectedStations)
    {
      case 0:
        // skip rejected station -> eqiv. to SignalStationsRange
        return !station.HasSignal();

      case 1:
        // confirmed noise station - problem: double rejections not covered
        return !(station.GetRejectedReason() == revt::eLargeTimeResidual ||
                 station.GetRejectedReason() == revt::ePolarizationDeselected ||
                 station.GetRejectedReason() == revt::eLonelyStation ||
                 station.GetRejectedReason() == revt::ePulseShapeDeselected);

      case 2:
        // all rejected stations will be treated as signal stations -> nothing to do
        return true;

      default:
        throw InvalidConfigurationException("invalid config fIgnoreRejectedStations.");
    }
  }

}